Five drugs were selected to be tested in combination (See Materials and Methods
). Each of the drugs reactivates KSHV with varying degrees. With the utilization of the five drugs that function in different yet potentially connected signaling processes (), KSHV reactivation can serve as an excellent model system illustrating how multiple cellular signals are processed. The five drugs are: Bortezomib, db-cAMP, Prostratin, Valproate, and Dexamethasone. Bortezomib is a proteasome inhibitor that at least in part reactivates KSHV by inhibiting NF-kB activity 
. DibutyrylcAMP (db-cAMP) is a cell-permeable cAMP analog that activates the PKA pathway 
. Prostratin activates the PKC pathway 
. Valproate shares structure and mechanism similarities with the histone deacetylase inhibitor butyrate 
. Dexamethasone is a glucocorticoid regulating the activation of some transcription factors and apoptosis-related genes 
Single-drug effects of KSHV reactivation and related cellular signaling.
In order to quantify the viral reactivation response, the RTA binding site in the PAN promoter was identified 
and a GFP reporter system was constructed. The reporter system BC-3-G, uses BC-3 cells (a primary effusion lymphoma cell line latently infected with KSHV) where a GFP protein is expressed under the control of a minimal lytic promoter of Polyadenylated Nuclear RNA (PAN), the most abundant KSHV early lytic transcript 
. Therefore, the expression of GFP following the activation of the PAN promoter served as a sensitive indicator of KSHV reactivation. The specificity of the reporter has been demonstrated in a previous study 
Measurement of virus reactivation was achieved using flow cytometry where we measured the number of activated cells, i.e., GFP positive, and the total number of cells, i.e., the number of dead and living cells. The reactivation rate (performance) of any given combination was set to be the ratio of GFP positive cells to the total number of cells including dead cells.
Modeling of mutli-signal induced KSHV reactivation
Investigation of the combinatorial effect of multiple participating pathways on reactivation can be achieved by treating the latently-infected cells with related chemical agents. Single drug dose curves for each chemical agent were obtained to determine the range of effectiveness of each individual chemical agent (). Based on the sensitive range of each individual agent determined from the curves, we selected the ranges of the concentrations to be used. Subsequently, the ranges were divided into ten concentrations using two-fold dilutions and setting the lowest concentration to zero (). The ten different concentrations of each drug comprised an input space of
possible drug combinations in total. Testing this number of combinations poses significant challenges (cost, labor, time, etc…). The choice of 10 concentrations depends on the shape and smoothness of the response and can be increased for finer sampling of the system response. However, the increase will lead to an increase in the number of tests.
Table of drug concentrations used in this study.
Using a uniform probability distribution over the set of all combinations of concentrations of five drugs, we randomly selected 600 different combinations to be experimentally tested. Six sets of experiments were conducted. In each set, 100 data points along with a positive control, a single drug (TPA) known to reactivate the virus 
, were evaluated using the GFP reporter system. The latently infected BC-3-G cells were treated with the combinations for one hour, after which the drugs were washed out, and measurements were taken 16 hours later to allow enough time for GFP synthesis and assembly upon reactivation.
The inputs (drug combinations) and their corresponding measured outputs (reactivation rates) were used to generate a mathematical model, KSHV reactivation model. The predictive reactivation model approximates the KSHV reactivation rate as induced by a combination of drugs within the specified range of drug concentrations. The model can be used to simulate and predict reactivation rates in response to all combinations of the five chemical agents. The combinations are not limited to the 600 tested combinations and include all combinations of the lower order mixtures, i.e., two, three, and four-drug combinations. The use of this relatively small number of combinations is facilitated by the assumption that the response function to these five drugs is reasonably smooth. If the response function is not very smooth then it would require testing of additional combinations to improve the accuracy and prediction power of the model.
Several methods can be used to generate a mathematical model. We utilize neural networks, linear regression 
, and partial least squares regression 
. Artificial neural networks are biologically inspired adaptive information processing systems. Artificial neural networks have been successfully applied to a wide range of problems in various disciplines including biological, medical, engineering, and financial 
. The combination of linear regression with partial least squares or all subset regression provided the ability to reduce the dimensionality of the problem and provides insight into which variables have the most influence on the observed responses.
We trained a multi-layered perceptron with the data and obtained a representative predictive model (see the Material and Methods
section). The model gave a correlation coefficient of more than 95% between the calculated and experimental data of the training set (). This indicates that the model has a reasonable prediction power but requires generalization. Therefore, the model was tested with an independently and randomly selected data set of 48 different combinations experimentally tested several months after the 600 points. The model was able to predict the corresponding reactivation rates with a correlation coefficient of 82%, a good fit considering the variability of cell responses due to varying cell conditions at different measurement times ().
Predictive modeling of reactivation rates.
Model-based optimization of KSHV reactivation
The predictive model generated provides the ability to determine combinations that can lead to high reactivation rates as predicted by the model. The simulated reactivation rates of all
combinations were enumerated. A simple sorting algorithm was used to rank the combinations in order of simulated reactivation rates. It is important to note that while a single best performing combination can be selected based on enumeration of all performances, the relevance of this best performing combination is not high due to measurement noise and modeling errors. Therefore, one is interested in looking at the distribution of top performing combinations. The top ranking 50 combinations were determined (). The distributions of individual concentrations within this group of points shows that lower to middle concentrations of Bortezomib are predominant. The distribution of concentrations of the other four drugs indicates that medium to high concentrations are predominant. The distribution of the performances within the top performing points indicates that the variation is within 3% of the maximum.
Characterization of the effect of drug combinations on KSHV reactivation.
An alternate approach to determine the top performing combinations is to utilize a search algorithm, deterministic or stochastic. Examples include gradient descent algorithms 
, genetic algorithms 
, the cross entropy method (CE) 
, as well as other stochastic search and combinatorial optimization algorithms. While a simple sorting algorithm suffices to sort all the performances, we apply a stochastic search algorithm here to search for optimal combinations based on the model to mimic similar experiments that we performed. This enables us to compare the outcomes of the two search experiments and to assess the possibility of running such algorithms to drive a set of experiments.
The cross entropy algorithm was implemented in silico using the KSHV predictive reactivation model (see Materials and Methods
). The simulated CE optimization showed that generally after about 14 iterations, the individual drug concentrations converged to 0 or 1.25 nM for Bortezomib, 4 mM or 8 mM for db-cAMP, 40 uM or 80 uM for Prostratin, 6 mM for Valproate, and 100 nM or 200 nM for Dexamethasone, to achieve consistently high reactivation (). The approach for optimizing combinations through the simple selection of the maximum possible dose of each drug does not result in a better reactivation rate than the optimized combination. The reactivation rate with the maximum doses of each agent, applied for one hour, is nearly 5%, where as it is about 42% with the optimized combination.
Experiment-based optimization of KSHV reactivation
An alternate approach to optimizing drug combinations is the use of a search algorithm implemented experimentally rather than on a mathematical model. Recently, several examples of this approach has emerged in biology 
. This approach can identify potent combinations and is useful in many situations where one is only interested in knowing which combination maximizes a pre-defined performance function.
The experimental cross entropy implementation proceeded in a sequence of experimental iterations. Our results showed that after 12 to 14 iterations, the drug concentrations converged to the ranges leading to consistently high reactivation rates (). The concentration ranges were 0–5 nM for Bortezomib, 4–8 mM for db-cAMP, 20–40 uM for Prostratin, 1.5–3 mM for Valproate, and a wide range of 0–200 nM (centered around 100 nM) for Dexamethasone. We further narrowed down the drug concentration ranges through another small set of iterations with drug concentrations more densely distributed within the initially determined ranges (). As expected, more consistently high reactivation rates were observed with the progress of the CE iterations. The optimal drug combinations obtained from the experimental CE method were consistent with the results from the simulated CE method as well as the direct enumeration (). This result experimentally validated the feasibility of a model-based approach in characterizing and optimizing multi-drug combinations.
Functional validation of selected combinations
Using two different approaches, we were able to identify a range of concentrations for which high virus reactivation rates are achievable. The results of the two approaches were consistent. To further validate the findings, we conducted sets of experiments to compare the performance of a selected combination from the identified range to the performances of single drugs.
The KSHV early lytic protein K8 is activated by, and expressed after the expression of KSHV RTA (ORF50). It is important for initiating viral DNA replication in the lytic cycle, thus a good marker for viral lytic replication. Western blot analysis of K8 showed that the selected drug combination can cause a much higher induction of K8 than any single drug. The conclusion was consistent 8 hours and 12 hours post treatment ().
Experimental validation of results.
Additionally, we looked at the effect of the selected drug combination on the KSHV lytic transcripts RTA (ORF50) and PAN. RTA plays a central role in regulating the switch from latency to lytic replication in KSHV 
. The activation of RTA (ORF50) is the first event in KSHV reactivation. It encodes the initiator of the viral lytic gene expression program. PAN (polyadenylated nuclear RNA), is the most abundant transcript made during the lytic cycle, and is directly induced by RTA 
. Quantitative analysis of these two lytic transcripts shows results similar to the western blot study of K8. Both lytic transcripts were induced approximately ten folds higher using the selected combination than the best concentration of any single drug. The results were also consistent at two different time points (). Our data shows that the combination treatment can potentially accelerate the reactivation process. Furthermore, we tested the virion production using Q-PCR upon treatment with a single drug and the optimal combination. The results show that there is an increase in virion production with the optimal combination over any single drug ().
Examining drug interactions
The signaling network involves complex connections between various molecules that can be perturbed through a large number of external signals. The signals can cause inhibition of certain molecules/pathways and stimulation of others. The interactions amongst these molecules or pathways are very complex and very hard to predict. Alternately, looking at interactions between the input signals and the measured cellular outputs can shed some light on the induced behaviors at the systems level. Particularly, we can uncover some of the interactions of the signaling that are involved in generating the responses upon stimulation with multiple stimuli.
Based on the predictive reactivation model, we simulated the interactions generated by the five drugs. The data represents a complex multi-dimensional data set. While some mathematical tools can be useful in reducing this complexity, one might be interested in visually examining the behaviors represented by such a large data set. To that end we created an interactive webpage which displays the KSHV reactivation rates for varying concentrations of the considered drugs 
Our findings indicate that the dose dependent effect of the individual drugs on reactivation greatly depended on the amounts of the other drugs within the same treatment (, webpage on accompanying CD). The results clearly indicate that drugs can interact to produce higher levels of cellular activity. However this improvement in reactivation is dependent on the concentrations of the drugs and needs to be optimized. The KSHV reactivation rate in the absence of drugs Valproate and Dexamethasone are less than the corresponding rates when these two drugs are present at certain concentrations (). The non-optimized addition of drugs to the system might not result in a noticeable improvement. In addition, the presence of appropriate doses of the drugs Valproate and Dexamethasone results in an increase of the effective range (the range for which high reactivation rates can be achieved) of drugs Bortezomib, db-cAMP, Prostratin. This provides the ability to use the drugs with lower concentrations while maintaining high reactivation rates.
Multi-drug response maps of KSHV reactivation.
The addition of low concentrations of Bortezomib to combinations of db-cAMP and Prostratin does not result in a significant increase in performance. Higher concentrations of Bortezomib result in a significant decrease in performance. Examining the effect of only adding Valproate to combinations of Bortezomib, db-cAMP, and Prostratin, we notice an increase in performance, indicating that Valproate interacts positively with the three-drug combinations to improve the reactivation. The sole addition of Dexamethasone to combinations of Bortezomib, db-cAMP, and Prostratin results in a smaller increase in performance. The increase becomes less when high concentrations of Bortezomib are used.
The above results reflect visual analysis of the responses, in the sequel, we seek to quantitatively analyze these interactions to determine the most significant ones. Such can be achieved using mathematical modeling similar to what is used for optimization. While a drawback of neural networks models is that they are black-box models and do not shed light onto how the different inputs are processed to produce the outputs, other modeling techniques can help in this regard. We fitted a linear regression model to represent the relationship between the drugs and the reactivation (see methods
section). The model utilizes 31 variables (regressors) that represent drug concentrations as well as interaction terms between the drugs.
The correlation coefficient between the experimental data and predicted data based on the linear model was 85%, the correlation coefficient for the additional 48 points was 83%. The model provides an insight into which factors play the biggest role in the response (). In agreement with the observations in the single dose-response curves and the neural network model, Prostratin and db-cAMP strongly influence virus reactivation. Additionally, there are other two and three-drug interactions that influence the response. Given this large number of model variables, we sought to find the key variables that affect the response. A partial least squares regression shows that around 10 components are sufficient to describe the variance in the output data (). The components of partial least squares model would be hard to interpret given the large number of variables. Instead, we pursue a subset selection algorithm based on all the possible subset regressions 
. The algorithm provides the best models of
variables. In total, the algorithm provides the best 31 models out of
Dimensionality reduction of the predictive reactivation model.
The residual sum of squares of the best models shows that there is no significant reduction in the residual sum of squares for models with more than 10 variables (95% reduction in the residual sum of squares). This indicates that 10 variables are sufficient to generate a model with comparable prediction and error to the 31-variable model (). The 10 variables of the best 10-variable model include the concentrations of the five drugs and products of two, and three-drug concentration (). This shows that the response is not only influenced by the individual drugs, but also by two and three-drug interactions. Most notably, there is strong negative interaction between Prostratin and Bortezomib, and strong positive interactions between db-cAMP and Prostratin, db-cAMP and Dexamethasone. A three-drug negative interaction between db-cAMP, Prostratin, and Dexamethasone is also present. Examination of models of 12 and 15 regressors shows that other three and four-drug interactions are present such as Bortezomib–Prostratin–Valproate, Bortezomib–db-cAMP–Prostratin, db-cAMP–Prostratin–Valproate, Bortezomib–db-cAMP–Prostratin–Valproate, and Bortezomib–db-cAMP–Prostratin–Dexamethasone.
Evaluation of effective subsets of combinations
Testing a system with five drugs provides advantages over studying mixtures of a smaller set of drugs. A system level study of combinations of multiple drugs enables fast and effective selection of a smaller subset of drugs that is most potent. Although a set of five drugs was used in this study, it is sometimes desirable to use a smaller number of drugs that can interact in a desirable way. We computed the maximum predicted reactivation rate for all possible mixtures of two, three, four, and five drugs, as well as for single drugs (). The figure shows there are significant differences between the maximum achievable reactivation rates using two, three and four drugs.
Evaluation of combinatorial effects of drugs on reactivation and cellular signaling.
For two-drug mixtures, there is over a six-fold difference between best and worst two-drug combinations. A mixture of Bortezomib and Dexamethasone or Valproate and Dexamethasone perform poorly even compared to a single drug. In contrast, a combination of db-cAMP and Prostratin have a reactivation rate higher than the sum of the individual reactivation rates. Prostratin and Valproate exhibit a similar behavior. This is consistent with our findings of strong positive interaction between Prostratin and Valproate. A mixture of Bortezomib and Prostratin does not improve on the best reactivation rate of Prostratin, suggesting negative or no interaction between the two drugs. This is also consistent with the findings presented above.
For three-drug combinations, the best combination is more than twice as effective as the worst mixtures. Furthermore, the best three-drug mixture is about 130% more effective than the best two-drug mixture. Four and five-drug mixture are slightly more effective than the best three-drug mixture. The four-drug mixtures generally perform better than than the three-drug mixtures.
Without a study of the combinations of five drugs, evaluating the reactivation rates for combinations of two drugs requires conducting 10 experiments individually to determine the maximum reactivation rate of the 10 possible two-drug combinations out of a set of possible five drugs. Selection of three or four-drug combinations requires similar experiments. Therefore, the combinations of the systems approach, computational tools, and experimental design enabled efficient multi-signal control of cellular/viral processes.