Dynamics of the Gene Regulatory Network of HIV-1 and the Role of Viral Non-coding RNAs on Latency Reversion

The use of latency reversing agents (LRAs) is currently a promising approach to eliminate latent reservoirs of HIV-1. However, this strategy has not been successful in vivo. It has been proposed that cellular post-transcriptional mechanisms are implicated in the underperformance of LRAs, but it is not clear whether proviral regulatory elements like viral non-coding RNAs (vncRNAs) are also implicated. In order to visualize the complexity of the HIV-1 gene expression, we used experimental data to construct a gene regulatory network (GRN) of latent proviruses in resting CD4+ T cells. We then analyzed the dynamics of this GRN using Boolean and continuous mathematical models. Our simulations predict that vncRNAs are able to counteract the activity of LRAs, which may explain the failure of these compounds to reactivate latent reservoirs of HIV-1. Moreover, our results also predict that using inhibitors of histone methyltransferases, such as chaetocin, together with releasers of the positive transcription elongation factor (P-TEFb), like JQ1, may increase proviral reactivation despite self-repressive effects of vncRNAs.


INTRODUCTION
Combined antiretroviral therapy (cART) is currently the most effective approach to control the chronic infection of HIV-1. However, cART does not eliminate the virus even with treatment intensification (Dinoso et al., 2009). This occurs because HIV-1 is able to form long-lived reservoirs by remaining latent within resting memory CD4+ T-cells (Siliciano et al., 2003;Siliciano and Greene, 2011;Cohn et al., 2015). Recently it has been proposed the use of LRAs in combination with cART to eliminate latently infected cells. Ideally this "shock-and-kill" strategy could purge viral reservoirs because when LRAs reactivate latently infected cells, those cells may be eliminated by self HIV-1 replication or by action of the immune system while cART prevents the formation of new viral reservoirs (Deeks, 2012). Despite many in vitro observations suggest that this strategy can be a promising approach (Deeks, 2012), clinical trials with LRAs have shown that it is ineffective in vivo (Bullen et al., 2014). Stochastic modeling of latently infected cells indicated that the clinical underperformance of LRAs is due to their inability to minimize the size of the viral reservoirs (Hill et al., 2014). Furthermore, this study suggested that LRAs must reduce the size of viral reservoirs 10,000-fold to prevent HIV rebounds after cART (Hill et al., 2014), an objective that cannot be reached with current treatments (Cillo et al., 2014).
The "shock-and-kill" strategy is based on the assumption that proviral reactivation depends only on the immunological activation of the infected cells. However, recent findings suggest that this assumption is not entirely true, since it has been observed that the provirus is able to autonomously regulate its latency using the positive feedback loop of transactivator of transcription (Tat) independently of cell activation (Razooky et al., 2015). During early stages of infection, Tat is synthesized at low levels that fluctuate because of cell's downregulation of the provirus (Weinberger and Shenk, 2007). When these transcriptional fluctuations are sustained, the activity of Tat initiates a positive feedback loop which boosts proviral transcription by recruiting P-TEFb in order to increase the synthesis of full-length viral RNAs (Weinberger and Shenk, 2007;Romani et al., 2010). In a biological context the two classical functions of positive feedback loops are to amplify and to sustain gene expression (Zhang Q. et al., 2014), however the architecture of the Tat circuit only amplifies transcriptional fluctuations making the gene expression of provirus transitory (Weinberger et al., 2005;Weinberger and Shenk, 2007). This architecture constitutes a mechanism of negative self-regulation of HIV-1, which may hinder viral reactivation (Razooky et al., 2015), and therefore may obstruct the activity of LRAs. Nevertheless, Tat is not the only structural component of HIV-1 that has a regulatory circuit. It has been observed that several vncRNAs have their own positive and negative feedback loops that may increase or suppress gene expression of the provirus (Groen and Morris, 2013;Saayman et al., 2014;Zhang Y. et al., 2014;Suzuki et al., 2015). It has been suggested that those vncRNAs have a secondary role on latency maintaining (Suzuki et al., 2015) and it is not clear whether such viral components participate in the low efficiency of the LRAs.
Current mathematical models of HIV-1 biology have been focused on transmission dynamics, posttreatment control, Vorinostat, and Romidepsin treatments, as well as the relation between reservoir size and reactivation (Hernandez-Vargas, 2017). However, none of these models addressed whether exist other paths to manipulate molecular components of the HIV to enhance latency reversion. Here we used Boolean and ordinary differential equations (ODEs) models to analyze the dynamics of the GRN of provirus to investigate how to reactivate more efficiently viral reservoirs with LRAs. In this network we included the interactions mediated by early viral proteins, vncRNAs, and epigenetic factors that regulate latency in resting CD4+ Tcells (Figure 1). It is important to remark that we used two different mathematical models in order to obtain results that represent the real dynamics of the GRN, independently of the model type chosen. The discrete model was used to calculate global properties of the network (attractors and its basins). The continuous model was used to measure changes in RNAs and protein expression levels of the GRN components. Both models consistently showed that the architecture of the GRN of wild type proviruses favors latency over activation state because of redundant interactions of vncRNAs. Furthermore, the models showed that reactivating effects of LRAs also stimulate the increase of vncRNAs, which reduces proviral protein expression. Finally, the models showed that the use of inhibitors of histone methyltransferases (HMTs) with releasers of P-TEFb, like chaetocin and JQ1 respectively, may increase proviral reactivation even in presence of vncRNAs.

MATERIALS AND METHODS
This work was performed in four stages: (1) Defining the GRN and its models, (2) Mathematical analysis of the GRN models, (3) Perturbation analysis of the models, and (4) Validation. The complete flux diagram of the methodology of this work is shown in Figure 2. During the first stage we constructed the GRN as well as the Boolean and the continuous models, then both models were analyzed separately. For the Boolean model, it was calculated its attractors with their respective attraction basins, then it was calculated the activation trajectory of the GRN and finally, it was evaluated the sensitivity of the model with the Derrida Test. On the other hand, it was calculated the equilibria of the ODEs model and it was evaluated the behavior of trajectories around such points with the analysis of stability, it was then evaluated the effect of particular changes in parameters values with the bifurcation analysis and finally it was evaluated the sensitivity of the model with a global sensitivity analysis. In the third stage it was performed a screening assay to find perturbations that reactivate latent proviruses and it was analyzed the dynamical features of such perturbations with discrete and continuous models. Finally, we validated both models with experimental data available from literature. In the following paragraphs of this section we present details of the protocols used in this work.

Construction of the Network
The GRN was built by compiling information from the literature on the molecular mechanisms that regulate HIV-1 latency inside resting CD4+ T-cells (Figure 1). This GRN included the main interactions of antisense long-non coding RNAs (asRNA), viral small interfering RNAs (vsiRNA), viral small activator RNA (vsaRNA), Tat, Rev, Nef, Vpr, and cellular factors that control gene expression of latent proviruses such as histone deacetylases (HDACs), histone acetyltransferases (HATs), and HMTs.
We incorporated to the GRN the most important molecules and viral components involved in the regulation of provirus gene expression, namely: the concentration of NF-κB, HATs, and HMTs; the activity of viral promoters 5 ′ LTR and FIGURE 1 | Gene regulatory network of HIV-1 provirus. Inside latently infected resting CD4+-T cells, the provirus can be activated by transcription factors and chromatin remodeling machinery of the host during immunological stimulation. When this occurs, the provirus expresses a wide variety of RNAs such as spliced viral mRNAs, as well as non-coding viral RNAs such as vsiRNAs (viral small interfering RNAs) and asRNA (long anti-sense RNA). The non-coding RNAs together with early proteins Tat, Rev, Nef, and Vpr have regulatory functions of the provirus either inducing or repressing viral transcription. Once the intracellular conditions are favorable to viral proliferation, late proteins like gp41, p24Gag, and other structural viral proteins are produced. In this figure the repressive interaction of this network are represented by pink T-bars, and activating interactions are represented by black arrows. Early and late proteins are shown in light and dark blue, respectively. Components of the host's transcription machinery are shown in purple, yellow, and black.

DISCRETE MODEL
For the discrete dynamics, the state of the nodes of the network in Figure 1 are represented by a set of binary variables, = {σ 1 , · · · , σ N }, each one taking the value 1 for activation and 0 for inactivation. The value of each variable σ n is determined by its k n regulators, denoted by {σ n 1 , · · · , σ n kn }, through the equation where f n is a Boolean function that depends on k n arguments ( Table 1). This function is constructed according to the inhibitory or activating nature of the interactions between σ n and its regulators (Kauffman, 1969). The discrete time t advances in integer steps; the time t ′ at which the state of the regulators is evaluated is such that t ≤ t ′ < t + t, where t is the time it takes to σ n to respond to a change in its regulators. Traditionally, Equation (1) is implemented simultaneously (synchronously) on all the nodes of the network. In this synchronous case t ′ = t and t = 1. In addition to the synchronous update, we also implemented two other updating schemes: asynchronous and semi-synchronous.
In the asynchronous scheme a permutation with repetition of the network nodes {σ 1 , · · · , σ N } is chosen. Let us denote as P = {σ p 1 , σ p 2 · · · σ p L } this permutation, where L ≥ N. Then at each time step t the nodes of the network are updated one by one following the order of this permutation: first σ p 1 at time t ′ = t+ 1 L , then σ p 2 at time t ′ = t + 2 L , and so on until σ p L is updated at time t ′ = t+1. When σ p i is being updated, Equation (1) is applied with t = i L and t ′ = t + i−1 L . After all the nodes in the permutation have been updated, the time t advances one unit and the process is repeated until an attractor is reached.
For the semi-synchronous scheme the set of all network nodes = {σ 1 , · · · , σ N } is partitioned into S subsets {M 1 , · · · , M s } such that S j=1 M j = . All the nodes contained in M j are updated synchronously, but the subsets {M 1 , · · · , M s } are updated asynchronously: the nodes in M 1 are updated at time t ′ = t + 1 S , the nodes in M 2 are updated at time t ′ = t + 2 S , and son on until the nodes in M S are updated at time t ′ = t + 1. When the nodes in M i are being updated, Equation (1) is applied with t = i S and t ′ = t + i−1 S . A full time step to go from t to t + 1 consists in the updating of all the subsets {M 1 , · · · , M s }, one by one in successive order. The construction of the permutation P for the asynchronous scheme and the subsets {M 1 , · · · , M s } for the semi-synchronous one was based FIGURE 2 | Flux diagram of methodology. This work was executed in four stages which are: (1) Defining the GRN and its models, (2) Mathematical analysis of the models, (3) Perturbation analysis of the models, and (4) Validation. In the figure line blue denoted the first stage, line purple corresponds to stage 2, line red is for stage 3, and finally gray line is assigned to the fourth stage. Pink color is used to represent the methodological steps made with the discrete model and blue color for continuous model. on biological phenomenology that reflects the way in which the activation cascade across the network may occur, and it is presented in the Supplementary Material.
It is well-known that the size of the basin of attraction is modified by updating scheme (Gershenson, 2002). The belonging of a network state to a particular basin of attraction strongly depends to updating scheme chosen. This has a biological equivalence, because the cellular environment is noisy and the order of gene expression may occur in different ways. However, there are some network states that always belong to same basin of attraction independently of updating scheme used. We call to this property as robustness under updating scheme. We hypothesize that the set of network states with this property are relevant for the biological behavior of the provirus. We call this set of states as intersection of the network states. We calculated the intersection of the synchronous, semi-synchronous, and asynchronous to determine the trajectory of activation of the provirus.

Stability of the Boolean Model: Derrida Map Test
The discrete model can exhibit two dynamical regimes, ordered and chaotic, and a phase transition between them, the socalled critical point (Aldana, 2003). The characterization of these regimes is given by the behavior of the avalanche of perturbations (produced by stochastic fluctuations, gene knockout, or gene over expression). In the chaotic regime, small perturbations spread throughout the network over time, producing big changes in the network state. Therefore, a network operating in a chaotic regime and submerged in a noisy cellular environment would have very unstable phenotypes. In the order regime, the perturbations die out over time, preventing the network to respond to new changing environmental conditions. In the critical point, the perturbations neither spread to the entire network nor disappear. They typically remain confined within a small fraction of genes. In order to characterize the dynamical regime, we define the normalized Hamming distance h(t) at time t between two network states as: In this equation σ n (t) is the state of the nth gene at time t in a trajectory starting out from a given initial condition, andσ n (t) is the state of the same gene in a different trajectory generated from a different initial condition. The Hamming distance h(t) can be considered as the normalized size of the avalanche of perturbations generated by differences the two initial conditions. The Derrida map h(t + 1) = M(h(t)) (Derrida and Pomeau, 1986) relates the size of the avalanche at two consecutive time steps. It can be shown that M(h) is a monotonic increasing function with the property that M(0) = 0 (if there is no perturbation at time t, there is no perturbation either at time (t + 1). The slope S at the origin of M(h) is the parameter that characterizes the asymptotic value of the Hamming distance, and hence the network dynamics. S is called the average network sensitivity. When S < 1 the network is operating in the ordered regime. If S > 1, the network exhibits chaotic behavior. If S = 1, the network is at the critical point. An intuitive definition (Krawitz and Shmulevich, 2007) is that S is the average fraction of genes that change their state at time t + 1 when a single gene is perturbed at time t (Supplementary Material). Therefore, to determine the stability of the network dynamics under perturbations in the initial conditions, one has to compute the network sensitivity S from the Derrida map M(h). Additionally, one can compute the network stability under permanent perturbations. We implemented two types of permanent perturbations: inhibition and overstimulation. For this, we set the state of one node, say σ j , equal to 0 or 1 all the time (regardless of the state of its regulators). Setting σ j = 0 for all time is equivalent to permanently inhibit this node, while setting σ j = 1 all the time is equivalent to having this node being constantly expressed. Let us denote as S j the network sensitivity when σ j is permanently perturbed (either inhibited or overstimulated), and as S 0 the sensitivity of the wildtype network. In order to compare the dynamical properties of perturbed proviruses vs. the WT provirus, we define the difference of sensitivity S as: This quantity measures how the network dynamics changes when one of the nodes is permanently perturbed. We performed the same type of analysis for the case in which two nodes σ i and σ j are simultaneously perturbed in a permanent way, either inhibiting or overstimulating them. This allows us to determine whether between-node epistasis exists that can modify the dynamics of the GRN.

Probability of Viral Activation
It is important to note that in the three updating schemes presented here, i.e., synchronous, asynchronous, and semisynchronous, the network dynamics are deterministic (both the permutation P and the subsets {M 1 , · · · , M s } are fixed). Therefore, in any of these updating schemes, after a transient time the network will fall into an attractor (a periodic pattern of activity). Several attractors may exist, and all the initial conditions that eventually fall into the same attractor are known as the basin of attraction of that attractor. As we show in the Results section, the HIV-1 network has several attractors. In some of them the network dynamics correspond to an active virus (the viral proteins are expressed, particularly p24Gag), whereas in the other attractors the dynamics correspond to an inactive virus (i.e., in the latency state with no expression of p24Gag). We refer to the former as the active attractors and to the latter as the inactive attractors. In order to determine the probability that a given initial condition leads to the active viral state, we compute the relative size of the activation state (W on ) by adding the size of the basins of attraction for all active attractors and dividing this sum by the total number of network states: where = 2 N is the total amount of network states, and |B(a k )| is the size of the basin of attraction of the k-th active attractor. Similarly, the relative size of latency state (W off ) was calculated as follows: These metrics determine the frequency of each state of the GRN that leads to an active or inactive attractor.

CONTINUOUS MODEL
In the continuous model, we represent the state of the nodes of the network in Figure 1 by the continuous variables {x 1 , · · · , x N }, which satisfy the general equation of mass balance ( Table 2) where the sums k J i n k and j J o n j represent all the fluxes that contribute to increase and decrease x n , respectively. The fluxes are presented in detail in Table 2, and the kinetic parameters (which were obtained from the literature), in the Supplementary Material. The Runge-Kutta 4-5 method was used to solve the system of ODEs. For extracellular pulses of TNF: For signals that stimulate HATs activity: For signals that stimulate HMTs activity: The discrete and continuous models present activation and latency states for WT proviruses and deletions of nef and vpr but only present attractors and equilibrium points for latency state when tat, p5 ′ LTR and splicing sites. This behavior is the same as reported for defective p ′ 5LTR mutants (Ho et al., 2013) and the splicing sites (Purcell and Martin, 1993), deleted tat (Verhoef and Berkhout, 1999), vpr (Rücker et al., 2004), rev and nef (Churchill et al., 2007). (B) Validation of the Boolean model. In this panel is shown the size of activation state (W on ) calculated with the synchronous, semi-synchronous, and asynchronous update scheme. In pink is shown increases of W on with respect to WT provirus. In the column of in vitro observations, "+" represents that there was an increase of viral reactivation because of the treatment and "0" indicates that there were no changes. The data for HDACis was obtained from (Cillo et al., 2014), for P-TEFb releasers from (Li et al., 2013), the use of Antagomirs from (Zhang Y. et al., 2014), combinations of Bryostatin with P-TEFb releasers from (Laird et al., 2015) and combinations of Bryostatin with HMTis from (Bouchat et al., 2012). (C) Validation of the ODEs model. In this panel is presented the normalized data of unspliced viral mRNAs levels obtained with the ODEs model and the corresponding values obtained from patients treated with bryostatin (Bullen et al., 2014), panobinostat (Laird et al., 2015), and JQ1 (Laird et al., 2015) vs. their corresponding simulation. Pearson correlation between both data sets showed a positive linear relationship, p = 0.0291, r (3) = 0.9708, which supports the validity of the model. The standard error of linear regression was 0.1613.
In these equations T 1 , T 2 , and T 3 are the activation intervals of the input signals (Supplementary Material).

Stability Analysis
The stability analysis of the continuous system was performed using the indirect method of Lyapunov (Supplementary Material). This method starts solving the ODEs in order to find the equilibrium points of the system. Then the ODEs are linearized using the Jacobian matrix to calculate the eigenvalues for all equilibrium points (Supplementary Material). Positive eigenvalues correspond to unstable directions in the phase space, whereas negative eigenvalues correspond to stable directions. If all the eigenvalues corresponding to one equilibrium point are negative, then that point is stable.

Bifurcation Analysis
The bifurcation analysis of the ODEs model was performed by changing one by one the parameters of the model. We focused our attention on the dissociation constants of NF-κB, association and dissociation constants of viral proteins, and degradation constants of RNA's and viral proteins (Supplementary Material).
Then, each parameter was varied three orders of magnitude, up and down of their reference value and after that; MATLAB was used to calculate the equilibrium points of the system with their corresponding stability.

Global Sensitivity Analysis
The sensitivity of the model against random perturbations was evaluated by assigning a uniform distribution to each parameter in which their reference value was taken as the mean and the standard deviation was assumed to be 10% (Supplementary Material). Then, each distribution was randomly sampled to obtain a set of parameters that were used as the inputs to solve the equations of the model during 1,500 units of time. After 10,000 iterations of this process, the concentration of p24Gag was used as the system's output to analyze the behavior of the model in response to random parameter variation. In all the simulations we set TNF(t) = 0, I HATs (t) = 0 and I HMTs (t) = 0.

Simulating Mutants and Treatments
The behavior of mutant proviruses during the condensation of viral nucleosomes and T-cells activation was modeled by  reducing 10-fold the splicing rate of nuclear mRNA of 4 kb (s 1 ). The nucleosomal condensation was modeled by providing square pulses of I HMTs , and T-cell activation was modeled by increasing the value of the NF-κB activity rate (k 1 ). The temporal effects of treatments with histone deacetylase inhibitors (HDACis), PKC agonists, P-TEFb releasers, histone methyltransferase inhibitors (HMTis), and antagonist micro-RNAs (antagomirs) on the GRN dynamics were simulated as follows: to simulate the rise on acetylation due to HDACis, we increased two-fold the reference value of the parameters of HATs activity (k 5 ). The increase the NF-κB levels due to PKC agonists (Mehla et al., 2010), was modeled by increasing the value of NF-κB levels (k 1 ) of the ODEs system. Considering that P-TEFb releasers, such as the compound JQ1, enhance the function of Tat to sequester P-TEFb and activate provirus (Li et al., 2013), we modeled this type of LRA by increasing the parameter associated to Tat activity (α 5 ). The effects of HMTis and antagomirs were modeled by reducing two-fold the reference value of the parameters of HMTs activity (k 8 ), synthesis of vsiRNA (α 2 ), and asRNA (α 4 ).
Mutant proviruses treated with HDACis were simulated by a two-fold increase in the value of the parameter of HATs activity (k 5 ) as a pharmacological overstimulation and setting to zero the values of the parameters for synthesis of Tat (α 5 and α 6 ), Nef (α 8 ), and Vpr (α 9 and α 10 ) as gene knockouts. The inhibition of vncRNAs was simulated by reducing 0-, 2-, 20-, and 200-fold the value of the parameters for the synthesis of vsiRNA (α 2 ) and asRNA (α 4 ). All parameters cited in this paragraph are listed in Supplementary Material.
Analogously to Equation (3), we define E 0 and E j as the normalized concentration of p24Gag mRNA for the wildtype network and when σ j is perturbed, respectively. The difference is a measure of the effect on the viral activation of perturbing the node σ j in response to pharmacological treatments.

Validation of the Models
The discrete and continuous models compatibility to reproduce the behavior of HIV-1 GRN was qualitatively evaluated by comparing the dynamical states of each model to the in vitro dynamics of provirus genic expression. To perform this, it was calculated the attractors of the discrete model and the equilibrium points of the continuous model for the wild type GRN and mutated networks p5 ′ LTR (t) = 0, Tat (t) = 0, Vpr (t) = 0, Rev (t) = 0, Nef (t) = 0, and mRNA4kbN (t) = 0. Then, the attractors and the equilibrium points were classified in activation state or latency state according to their p24Gag expression level (i.e., latency state was assigned to attractors and equilibrium points that do not express p24Gag and activation state was assigned when p24Gag is expressed). These results were compared to in vitro observations reported for the wild type provirus, defective p ′ 5LTR mutants (Ho et al., 2013), deleted tat (Verhoef and Berkhout, 1999), vpr (Rücker et al., 2004), rev and nef proviruses (Churchill et al., 2007), as well as deletions on the splicing sites (Purcell and Martin, 1993; Figure 3A). Once compatibility of the models was proved, the discrete model was qualitatively validated by comparing the size of activation state (W on ) of the nodes perturbations HATs (t) = 1, Tat (t) = 1, vsaRNA (t) = 1 and the combinations (NFκB (t) = 1, HATs (t) = 1), (NFκB (t) = 1, HMTs (t) = 0), against their in vitro equivalences, which are treatments with HDACis (Cillo et al., 2014), P-TEFb releasers (Li et al., 2013), the use of Antagomirs (Zhang Y. et al., 2014), combinations of Bryostatin with P-TEFb releasers (Laird et al., 2015) and combinations of Bryostatin with HMTis (Bouchat et al., 2012). In Figure 3B is shown the outcome of this comparison, which pointed out that the discrete model is able to predict at qualitative level changes occurred on latency reversion reported in vitro. The continuous model was validated by comparing the levels of genomic RNAs obtained in silico against ex vivo data (Laird et al., 2015). To perform this, it was normalized the levels of p24Gag obtained with the ODEs model for making a linear regression analysis and Pearson correlation with 5% of significance (α = 0.05; Figure 3C). These analysis showed that there is a significant positive relationship between ex vivo and in silico data sets, R 2 = 0.9426 with p < 0.05, which suggest that the ODEs model is able to predict variations over concentration levels of molecular components of the HIV-1 GRN.

FIGURE 4 | Proviral activation is repressed by vncRNAs and rescued by Tat. (A)
Relative weight of the activation state (in percentage). In this panel is shown the result of calculating the probability of activation using the synchronous, semi-synchronous, and asynchronous update schemes. In all cases it is shown that latency is favored over activation once provirus is integrated in the host genome. (B) Trajectory of activation obtained from intersection of the three update schemes. This panel shows that the presence of NF-kB, p5 ′ LTR, Tat as well as the absence of vncRNAs (i.e., vsiRNA and asRNA) and HMTs are obligatory conditions to activate latent proviruses. (C) Sensitivity of the network calculated with Derrida's mapping test. The data obtained from the Boolean model suggested that Tat and vncRNAs are the main proviral regulators of latency and activation.

T-Cell Activation May Not Induce Expression of the Provirus
Razooky and coworkers found evidence suggesting that proviral latency is mainly regulated by the transactivation of 5 ′ LTR mediated by Tat instead of T-cell activation, which implies that latency regulation may be an autonomous process (Razooky et al., 2015). It is in the light of this finding that, the role of epigenetic factors on the performance of Tat's autonomous behavior was investigated. To accomplish this, we analyzed the attractors of the Boolean and its basins of attraction in presence of cellular signals that stimulate epigenetic regulators such as HMTs and HATs, and activators of the NF-κB pathway like TNF. In the three update schemes it was found 12 punctual attractors (the same in the three schemes) which were classified in two groups according to expression of viral proteins as follows: (1) attractors that produce late proteins like p24Gag (activation attractors); and (2) attractors that lack protein expression (latency attractors; see Table 3).
The Boolean model shows that the activation attractors can be reached with or without cellular stimulation of HATs and TNF (Table 3), which agrees with previous observations that demonstrate the persistence of provirus expression in resting CD4+ T-cells (Razooky et al., 2015). However, this dynamics always requires the absence of the silencing produced by the HMTs activity ( Table 3). The probability with which the provirus reaches latency and activation was investigated by calculating the relative size of the activation state (W on ) as well as the relative size of the latency state (W off ). It was found that W on is always smaller than W off (Figure 4A) even when transcription stimulatory signals like HATs and the NF-κB pathway are turned on. These results suggest that even in the context of T-cell activation, provirus may remain latent because of its autonomous dynamics, which is limited by epigenetic silencing.

Viral Non-coding RNAs Are Essential to Regulate Latency
Previous reports showed the importance of Tat as the unique virus-encoded regulator of HIV-1 autonomous behavior (Weinberger et al., 2005;Razooky et al., 2015). However, a virus-encoded siRNA that also promotes provirus activation has been found recently (Zhang Y. et al., 2014). Additionally, other virus-encoded regulators, such as vncRNAs that directly repress provirus gene expression have been found (Groen and Morris, 2013;Saayman et al., 2014;Suzuki et al., 2015). Therefore, the role of vncRNAs on the regulation of proviral latency was investigated by searching for common states in all basins of attraction of activator attractors obtained with the three updating schemes ( Figure 4B). Using this procedure we found a set of GRN states that abrogate latency ( Figure 4B). This set of states indicates a general pattern that results in provirus activation, which agrees with previous reports and requires: high levels of NF-κB (Westendorp et al., 1995), no epigenetic silencing by HMTs (Jordan et al., 2003;du Chéné et al., 2007), genomic integrity of provirus (Ho et al., 2013), high levels of Tat (Weinberger et al., 2005;Razooky et al., 2015), and the absence of repressive vncRNAs (denoted by asRNA and vsiRNA; Figure 4B). This result suggests that Tat and repressive vncRNAs are essential virus-encoded regulators of latency establishment and activation.

HIV-1 Is Resistant to Drugs and Intracellular Perturbations
Genetic networks of organisms are able to maintain and adapt their operation in response to environmental changes. Previous studies have shown that the coexistence of robustness and adaptability observed in genetic networks is characteristic of systems operating at the critical point, i.e., at the border of chaos and order (Balleza et al., 2008). This dynamical feature has been reported for genetic networks of A. thaliana, D. melanogaster, S. cerevisiae, E. coli, B. subtilis (Balleza et al., 2008) as well as of mice macrophages (Nykter et al., 2008). It has been suggested that criticality is essential to ensure the evolution of any organism (Balleza et al., 2008). We investigated the presence of critical dynamics in the HIV-1 GRN. To do this, the effect of massive perturbations on the GRN was evaluated using the Derrida mapping test. When the network sensitivity S for the provirus GRN was computed, it was obtained S = 1.0031 which means that the network operates in a critical regime ( Figure 4C). Therefore, this network shows equilibrium between robustness and adaptability in resting CD4+ T-cells ( Figure 4B). This result suggests that the regulation of the expression of the HIV genome is robust against intracellular perturbations and it can be adapted in response to chronic perturbations, such as those produced during cART or treatments with LRAs. It should be noted that the HIV-1 network has constructed taking into account the activating and inhibitory interactions reported in the literature without considering criticality as a relevant criterion. The result showing that the dynamics of the HIV-1 GRN is so close to criticality is unexpected.

The Architecture of the HIV-1 GRN Allows Viral Rebounds and Persistence
Previous observations on the dynamics of Tat's positive feedback loop demonstrated that this circuit is able to amplify transcriptional fluctuations of provirus by itself, and its activity tends to decay toward a latency stable state (Weinberger et al., 2005). It has been proposed that delays on Tat's activity facilitate latency establishment (Weinberger et al., 2005), which could maintain proviral reservoirs during cART . However, it is unknown whether other viral components like Vpr, Nef, and vncRNAs modify the dynamics of Tat's circuit. In this direction, we extended previous findings by analyzing the provirus gene expression dynamics in the presence of Tat and other viral interactions that regulate proviral transcription, such as those mediated by vncRNAs and positive feedback loops of Nef and Vpr (Varin et al., 2003;Liu et al., 2014). To this end, it was used the continuous model to analyze temporal variations of the dynamics of the levels of provirus proteins and RNAs. It was performed the stability analysis of the ODEs model with a set of reference parameters (Supplementary Material), and found two equilibrium points that correspond to activation and latency states, i.e., the levels of p24Gag were zero for latency state and distinct to zero for activation state (Table 4). In this regard, the stability analysis showed that the activation state was stable and the latency state was unstable ( Figure 5A). Then, the sensitivity of the system against fluctuations was evaluated by performing a global sensitivity analysis finding that the dynamics of the system was robust against perturbations (Supplementary Figure 1), and the mean value of p24Gag during activation state was 11.2 with a variance of 7.2. This suggests that once activation state is reached, the provirus expression is resistant to variations of the intracellular environment. We searched for parameters that change stability of the equilibrium points of the GRN performing bifurcation analysis. Indeed, it was found a transcritical bifurcation (Figure 5B) on the value of NF-κB activation constant (k 1 ), parameters related to splicing of viral mRNAs (s 1 ), and the activity of the 5 ′ LTR promoter (k b ). Bifurcation analysis showed that latency is stabilized when the values of these parameters are decreased ( Figure 5C). This observation is congruent with in vitro reports of conditions that stabilize latency, such as low levels of NF-κB (Westendorp et al., 1995), deficient splicing sites (Purcell and Martin, 1993), and deletions on 5 ′ LTR promoter (Dar et al., 2014) (Table 1). Then, we investigated the possible function of this bifurcation in the context of intracellular infection of HIV-1. To implement this, it was compared the performance of WT Phase portrait of the system around the equilibrium point corresponding to latency (black cross), and the temporal performance of the system are shown. In the phase portrait, trajectories of the system are repelled to activation state, on the temporal plot of the system; p24Gag reaches an expression stable state. This simulation was made with our reference value for NF-κB availability (k 1 ). (B) Transcritical bifurcation on the GRN. We found that variations on parameters related to availability of NF-κB, the activity of 5 ′ LTR promoter and the splicing of viral mRNAs change the dynamical behavior of the system. The critical parameters to obtain this bifurcation are included in Supplementary Table 5. (C) Stabilization of latency. When we decreased 10-fold NF-κB availability, all trajectories in the phase portrait of the system converge to latency state (black cross), in the temporal plot this can be observed as a transient activation of protein expression that eventually decays. (D) Biological role of transcritical bifurcation. In the absence of this bifurcation, defective proviruses decrease their ability to relapse after a period of repression. This simulation was made by decreasing 10-fold the splicing rate of nuclear mRNA of 4 kb (s 1 ); gray bars indicate nucleosome compaction due to HMTs activity. (E) Molecular origin of transcritical bifurcation. Individually, positive feedback loops of Tat, Nef, and Vpr have a transient activity (as observed in panel C), however, transcritical bifurcation emerges when loops are combined. nef, vpr, and tat were simulated by setting to zero the synthesis parameters of Tat, Nef, and Vpr (Supplementary Material). Collectively, these data suggest that redundant activation of NF-kB mediated by Tat, Nef, and Vpr ensures proviral reactivation after a period of repression.
provirus vs. mutated provirus that have attenuated splicing rates (10-fold lower of the reference value for s 1 ) in presence or absence of epigenetic silencing (i.e., when HMTs are active). It was found that transcritical bifurcation allows viral rebounds of the WT provirus after cellular inhibition (Figure 5D), which suggests that persistence may be "hardwired" on the HIV-1 genome. This assay shows that 51% of the perturbations permanently silence provirus' expression; where "reactivation" refers to perturbations that suppress latency attractors, "no changes" refers to perturbations that allow the coexistence of latency and activation attractors, and "permanent silencing" refers to perturbations that abrogate activation attractors.
On the other hand, proviruses that lack transcritical bifurcation can be easily controlled by the host's HMTs ( Figure 5D). These results suggest that the transcritical bifurcation of the provirus GRN may provide two dynamical behaviors: (1) for repressive transcriptional environments, such as during cART, the provirus latency will be stabilized allowing reservoirs maintenance, and (2) for non-repressive transcriptional environments, the provirus favors a strong activation in order to ensure the production of viral progeny and to counteract the intracellular silencing mechanisms. These properties may explain the viral rebounds after cART and why HIV-1 cannot be silenced by host.
The Activating Core of the GRN Consists of the Positive Feedback Loops of Tat, Nef, and Vpr The molecular basis of the transcritical bifurcation was investigated comparing the activity of intact provirus vs. the activity of mutant proviruses. Mutant proviruses were simulated in the continuous model by setting to zero all parameters related to the synthesis of viral proteins Tat, Nef, and Vpr. It was observed that the Tat's positive feedback circuit always produces a stable branch on latency state, which in biological terms is a transient activation followed by latency stabilization dynamics as reported by Weinberger et al. (2005) (Figure 5E). However, combining Tat positive feedback with Vpr and Nef produces the transcritical bifurcation, in which latency can be destabilized ( Figure 5B). We also observed that in the absence of Tat the remaining positive feedback loops were able to temporarily perturb latency during stimulation, producing transitory gene activation, but their effect was negligible compared to that observed in the presence of Tat ( Figure 5E). Thus, the transcritical bifurcation is sustained by all the positive feedback loops of the viral proteins Tat, Vpr, and Nef ( Figure 5E). Considering that all the positive feedback loops of HIV-1 promote NF-κB activation (Figure 1), it is reasonable to think that the redundancy on NF-κB stimulation is the cause of the transcritical bifurcation and its amplifying properties.

Permanent Stabilization of Latency Occurs More Frequently Than Reactivation
Recently, it has been proposed that compounds that increase fluctuations of transcriptional basal levels may enhance the performance of LRAs (Dar et al., 2014). Such compounds indirectly target the 5 ′ LTR promoter, increasing its activity. We extended this result by searching for sensitive interactions that could increase proviral reactivation in the presence of LRAs. To this end, it was used the Boolean model to explore all possible perturbations of the provirus GRN by combining inhibition and stimulation of the GRN nodes using a screening assay (Figure 6). It was found that 51% of the perturbations eliminated activation attractors, which suggests that those perturbations are able to induce permanent silencing of the provirus (Figure 6D). On the other hand, it was found that only 28 of the 648 theoretical perturbations can be performed in vivo using current LRAs and antagomirs (Table 5). Remarkably, some of these perturbations have not been tested yet. These results suggest that it would be easier to induce the permanent silencing of HIV-1 proviruses rather than reactivating them ( Figure 6D).

Inhibition of HMTs and Stimulation of P-TEFb Increases Proviral Reactivation
We then characterized the dynamical properties of 28 promising perturbations produced with LRAs and antagomirs ( Table 5).
To do this, the dynamical performance of each perturbation was compared to the dynamics of the WT provirus. It was used the Boolean model to calculate the relative size of the activation state (W on ) and the difference of sensitivity ( S). Similarly, it was used the ODEs model to determine the difference of p24Gag expression ( E) for each perturbation. It was found that all reactivation perturbations increased W on , except HATs(+) (Figure 7A) which is equivalent to using HDACis (Table 5). Moreover, all reactivating perturbations decreased network sensitivity ( Figure 7B) and the ODEs model showed that all perturbations, except HATs(+), increased the expression of p24Gag ( Figure 7C). Remarkably, the discrete model showed that inhibition of HMTs and overstimulation of Tat, i.e., HMTs(-), Tat(+) precludes latency attractors, which means that provirus is always active ( Figure 7A). Analogously, the ODEs model showed that HMTs(-), Tat(+) increases E to the maximum ( Figure 7C). It is important to note that the pharmacological equivalence of HMTs(-), Tat(+) can be implemented with HMTis and P-TEFb releasers (Li et al., 2013;Table 5). In Table 5 are shown the pharmacological treatment equivalent for the other latency reversing perturbations.

The Performance of LRAs Is Hindered by vncRNAs
Recent reports showed that HDACis are not effective to reactivate latent proviruses (Bullen et al., 2014;Cillo et al., 2014). In agreement with these reports, the models showed that HDACis do not produce changes in the activation state ( Figure 7A) and do not increase p24Gag expression levels ( Figure 7C). However, it has been reported that HDACis increase transcription of provirus (Mohammadi et al., 2014). To explain the HDCAis underperformance, the existence of unknown posttranscriptional mechanisms that counteract protein synthesis have been proposed (Mohammadi et al., 2014). Furthermore, it has been reported that HDACis like SAHA (Vorinostat) may increase the levels of cellular non-coding RNAs (Lee et al., 2009). Taken together these observations suggest that HDACis increase provirus transcription as well as the levels of viral and cellular non-coding RNAs, which contributes to silencing protein expression of provirus. We explored this hypothesis by comparing W on , S, and E for each HDACis perturbation with and without vncRNAs (see section Methods). It was found that the suppression of vncRNAs enhances HDACis performance, increasing the values of W on (Figure 8A), S (Figure 8B), and the expression levels of p24Gag ( Figure 8C). These data suggest that HDACis may promote the synthesis of vncRNAs, which may explain why these LRAs increase provirus transcription but not protein expression (Figure 8D).

Inhibition of vncRNAs Is Not Sufficient to Stimulate Proviral Reactivation
The results just presented indicate that inhibiting vncRNAs could enhance the effect of LRAs ( Figure 8D). However, it is not clear whether vncRNAs inhibition can also stimulate the reactivation of mutant proviruses. Therefore, we used the ODEs model to address this question and compared the expression levels of p24Gag in defective provirus treated with HDACis at different intensities of vncRNAs inhibition. It was found that mutant proviruses that lack the Tat protein can be reactivated to a lesser extent than intact proviruses (Supplementary Figure 2). However, defective proviruses that lack two or more positive feedback loops cannot be reactivated, even with the inhibition of vncRNAs (Supplementary Figure 2). These results suggest that inhibition of vncRNAs cannot ensure the total reactivation of proviral reservoirs.

DISCUSSION AND CONCLUDING REMARKS
The long-lived latent reservoirs of HIV-1 are the main barrier to eradicate it. Several efforts to purge viral reservoirs have been performed using LRAs, unfortunately none of them were effective in vivo (Bullen et al., 2014). Until now it is not known the causes of the underperformance of LRAs. In this work, we analyzed in silico the functioning of the provirus' gene expression in order to investigate the ineffectiveness of LRAs. To this end, we constructed the GRN of provirus and modeled its dynamics using ODEs and logic rules. Both models predicted that vncRNAs are the main negative regulators of the gene expression of provirus and they are also implicated in the underperformance of LRAs. Finally, both models predicted that treatments with HMTis and P-TEFb releasers are the best way to maximize latency reversion.
Traditionally it has been thought that Tat is the only virusencoded regulator of the HIV latency. However, recent evidence shows that vncRNAs are also essential to control proviral latency. Saayman and colleagues characterized an HIV-encoded long anti-sense RNA which its inhibition triggers reactivation in latently infected cells (Saayman et al., 2014). Zapata and coworkers showed that this long anti-sense RNA is able to silence the gene expression of provirus by stimulating HMTs (Zapata et al., 2017). Thus, we investigated the role of vncRNAs on the dynamics of provirus' gene expression. The first dynamical particularity of the GRN was that the weight of the latency state (W off ) was higher than the weight of the activation state (W on ), regardless the cell's activation state ( Figure 4A). After analyzing the set of intracellular environments that activate the GRN (Figure 4B), we noted that activation requires the presence of Tat and the absence of vncRNAs. Additionally, the inhibition of vncRNAs increased the W on ( Figure 8A). Taken together these results indicate that vncRNAs are the main negative regulators of the provirus' genic expression.
The next question to address was how vncRNAs and Tat operate together to regulate latency. Previous reports demonstrate that Tat's positive feedback loop has a strong transient activation that eventually decays to a stable latency state (Weinberger et al., 2005;Weinberger and Shenk, 2007). The same behavior was observed on the Tat's circuit of the FIGURE 9 | Molecular mechanisms of self-regulation of proviral latency. According to our results, activation state is mainly produced when Tat concentration reaches high levels. On the other hand, the provirus induces its latency when Tat's concentration is not optimal and the levels of vncRNAs are high.
GRN (Figure 5E), as well as in other positive feedback loops mediated by Nef and Vpr ( Figure 5E). Interestingly, we found that a transcritical bifurcation appears when these circuits were combined (Figure 5B), and such a bifurcation allows gene expression rebounds after long periods of repression ( Figure 5D). It seems likely that the Tat's circuit is enhanced by Nef and Vpr in order to overcome the downregulation of vncRNAs and the host. However, an uncontrolled enhancement of the gene expression of provirus could have negative effects on the viral reservoirs. Rouzine and colleagues found that a high rate of proviral activation avoids the establishment of latent reservoirs, which decreases the prevalence of HIV-1 . They also observed that fluctuations on the transient activity of Tat, decreases the frequency of provirus' activation which stabilizes viral reservoirs . Expanding these observations, our results showed that in addition to Tat's fluctuations, vncRNAs also reduce the activation of provirus. Thus, vncRNAs together with Tat's transient activity may be responsible for the chronic stabilization of latency, condition required to maintain the viral reservoirs (Figure 9). Furthermore, we investigated the role of vncRNAs on the underperformance of LRAs. The screening assay (Figure 6) showed that 28 perturbations of the GRN can be implemented with LRAs and antagomirs (Table 5), being the combination of HMTis with P-TEFb releasers the most prominent of all. However, perturbations made with HDACis did not increase protein expression of provirus (Figure 7), as reported by Cillo et al. (2014). Mohammadi et al found that HDACis only increase provirus' transcription but did not affect protein expression (Mohammadi et al., 2014). They proposed that this occurs because of post-transcriptional mechanisms that hinder protein expression (Mohammadi et al., 2014). In this direction, our results predicted that the levels of vncRNAs increased in response to HDACis (Figure 8). Hence, it seems likely that treatments with HDACis stimulate proviral transcription as well as vncRNAs, which eventually avoids protein expression. This hypothesis may explain the underperformance of treatments with LRAs reported in vivo.
The final question to address was how to enhance the performance of LRAs. The screening assay showed 28 feasible treatments to disrupt latency by using micro-RNAs and current LRAs ( Table 5). In this direction the treatment that maximizes the probability to reactivate proviruses (given by the value of W on ) uses HMTis and P-TEFb releasers ( Figure 7A). The action mechanism of this treatment consists in increasing Tat's levels with P-TEFb releasers while the activity of HMTs is blocked, which is the main downstream target of vncRNAs (Zapata et al., 2017). Therefore, blocking molecular effectors of vncRNAs and enhancing Tat activity is the best way to increase viral reactivation. It is of our interest to test the effectiveness of the treatments proposed in Table 5 with ex vivo cultures obtained from HIV patients, in order to determine whether such treatments could be promising for therapeutic implementation.
Nevertheless, our results also showed an interesting scenario that has a distinct approach to control HIV-1. The screening assay showed that 51% of perturbations permanently silence the provirus genic expression (Figure 6D). It is noteworthy to say that the most of perturbations that permanently silence the provirus, inhibit nodes related to proviral transcription such as p5 ′ LTR and unspliced, spliced and partially spliced viral mRNAs (Figure 6A). This implicates that HIV-1 can be permanently controlled by the induction of hypermutation of its genome. A possible mechanism to implement this strategy can be achieved with APOBEC3G, which is the enzyme that naturally hypermutates HIV-1 as a part of intracellular antiviral response. In this context, APOBEC3G is inhibited by Vif in order to allow the progression of HIV-1 infection. However, recent findings suggest that drugs that stimulates ASK1 (apoptosis signalregulating kinase 1) also restore the APOBEC3G function even in presence of Vif (Miyakawa et al., 2015). Thus, an alternative path to control HIV-1 infection may employ APOBEC3G inducers in conjunction with cART.
Current treatments to reactivate latent proviruses may fail because HIV uses its vncRNAs as negative regulators to maintain latency. Some LRAs like HDACis could increase the levels of vncRNAs, consequently reducing their effectiveness to revert latently infected cells. Our results suggest that the best treatment to avoid the repressive effects of vncRNAs is to use an HMTis like chaetocin, together with P-TEFb enhancers. Treatment that could have potential for efficient reactivation of the HIV-1 provirus should be clinically tested.

AUTHOR CONTRIBUTIONS
AB and CT-S conceptualized, designed, and performed all computational experiments of this study. RG and JD contributed designing experiments, interpreting, and supervising this study.
AB and CT-S wrote the first draft of the manuscript. RG and JD reviewed drafts and approved final version of the manuscript.

FUNDING
This work was supported by CONACYT, and by the PRODEP funding program of the Universidad Autónoma del Estado de Morelos. AB was supported by CONACYT through the Ph.D. Scholarship number 27667. CT-S was supported by a PRODEP Postdoctoral grant 103.5/14/11795.