A Stochastic Binary Model for Regulation of Gene Expression to Investigate Treatment Effects Targeting RKIP

In this manuscript we use an exactly solvable stochastic binary model for regulation of gene expression to analyse the dynamics of response to a treatment aiming to modulate the number of transcripts of RKIP gene. We demonstrate the usefulness of our method simulating three treatment scenarios aiming to reestablish RKIP gene expression dynamics towards pre-cancerous state: i. to increase the promoter’s ON state duration; ii. to increase the mRNAs’ synthesis rate; iii. to increase both rates. We show that the pre-treatment kinetic rates of ON and OFF promoter switching speeds and mRNA synthesis and degradation will affect the heterogeneity and time for treatment response. Hence, we present a strategy for reducing drug dosage by simultaneously targeting multiple kinetic rates. That enables a reduction of treatment response time and heterogeneity which in principle diminishes the chances of emergence of resistance to treatment. This approach may be useful for inferring kinetic constants related to expression of antimetastatic genes or oncogenes and on the design of multi-drug therapeutic strategies targeting master regulatory genes.


Introduction
Despite recent advances in our understanding of cancer biology, the use of quantitative methods to integrate the plethora of generated data to design treatments targeting metastasis is still in its infancy [1,2]. The genotypic variability and intrinsic randomness of biochemical reactions governing epigenetics underpin the multiplicity of cancer cell phenotypes commonly termed as tumor heterogeneity [3][4][5][6][7][8]. Additionally, cellular processes are controlled by several networks of chemical reactions characterized by changeable topologies and functional redundancies that equip cells with adaptation and robustness capabilities [9][10][11][12]. And, the numerous characteristic timescales of the chemical processes taking place inside the cell adds one more layer of cumbersomeness [13,14]. Thus, the enhancement of therapeutic strategies requires the analysis and the "engineering" of the dynamics of treatment response in a complex system composed by several interacting components having a multiplicity of characteristic time scales and subjected to randomness. Finding building blocks having those features may provide useful insights on how to orchestrate the dynamics of a such a complex system [15]. For that, the exactly solvable stochastic model for transcription of a binary gene [16][17][18][19] is a good candidate to be used as a prototype for simulating enhanced treatment strategies.
To develop our theoretical approach for therapeutic design we start selecting the gene encoding the Raf kinase inhibitory protein (RKIP) as a model system because it plays a major role on regulation of the dynamics of multiple components of a cell [20]. Because cancer lethality is mainly caused by metastasis, the choice of RKIP is promising as its concentrations are typically reduced in metastatic cancers [21][22][23]. That negative correlation turns RKIP and RKIP-related gene signatures into useful biomarkers of metastatic risk in cancer patients [24]. The characterization of RKIP as an anti-metastatic gene is reinforced by experiments demonstrating that its overexpression blocks in vivo invasion and metastatic progression [25,26]. Hence, we focus on possible strategies to increase the amounts of RKIP in cancer cells by re-modulating its expression profiles towards the pre-cancer regimens.
For that goal, we simulate the effect of application of multiple drugs each targeting a different kinetic rate of a two state stochastic model for gene transcription. Because of treatment, we will consider the kinetic rates of the model to be time-dependent in response to drug application. That enables one to estimate the speed and heterogeneity of response to treatment by, respectively, computing the dynamics of the average number of transcripts and their variance. The exact solutions of the model at constant kinetic rates have two characteristic time scales, one related with the promoter state switching and the second being the mRNA lifetime. Recently, we demonstrated that the ratio between those two time scales may be used to classify the qualitatively distinguishable noise regimens of gene transcription on the binary model [19,27] and the reliability of information about the promoter state that is transmitted by gene products [28]. Here, we assume that the addition of a drug cause the kinetic rates of the model to become time dependent and that this effect decays exponentially such that has up to four additional time scales for the problem. Then, we show that treatment targeting a gene having a fast switching promoter, and expressing as a quasi-Poissonian process, enables the fastest and least heterogeneous response to treatment. If the gene has a slow switching promoter, the response will be slower and have higher heterogeneity. A gene expressing in bursty fashion will enable a fast response having maximal heterogeneity. Then, we build upon the previous analysis to design an enhanced treatment enabling a faster response with reduced heterogeneity that is independent of the pre-treatment time scales. The enhanced treatment is based on reduced drug dosages to reduce the chances of toxicity and emergence of resistance caused by compensatory effects [1,2].
The remaining of this manuscript is organized as follows. The assumptions underlying the qualitative and quantitative models that we propose are presented on Section 2. The obtained results are shown on Section 3 and we discuss them on Section 4. Our concluding remarks are presented on Section 5.

Methods and Models
Gene expression is the key mechanism keeping cell function by means of an extensive molecular machinery which transforms the genetic information into molecular function. Gene expression can be described as a two stages process, namely, the gene transcription that generates the mRNA, and the mRNA translation which has proteins as their products. Here we focus on gene transcription and provide a simplified picture of how it happens in eukaryotic cells: the RNA Polymerase protein complex (RNAPolII) binds to a specific region of a gene, the promoter site, and begins transcription elongation. The binding of RNAPolII to the promoter site may be regulated by its interaction with a DNA region called enhancer. The enhancers are regions of the DNA that interact with the transcription factors that can provide a positive or negative regulation of the binding of RNAPolII to the promoter site. This is the process that we will use to model RKIP transcription.

A brief description of the molecular role of RKIP
RKIP protein is a regulator of kinases that directly binds to Raf kinase [29] and is involved on regulation of signalling pathways such as the Raf-Mek-Erk cascade and the NF-κB-related pathways [30,31]. Both pathways modulate cell proliferation, the former regulates differentiation, while the latter is responsible for inflammation, and anti-apoptosis processes. In addition, RKIP modulates other fundamental cell signaling pathways involving heterotrimeric G-proteins, keap1/nrf2, STAT3 and GSK [32]. Because of its interaction with multiple pathways, we consider RKIP as a master modulator of cellular processes. The amounts of RKIP proteins inside the cell can be regulated at multiple steps of expression, from pre-transcription of the gene to post-translation [33]. A plethora of metastatic solid tumors have RKIP downregulated or lost and experimental data suggests that that happens because of transcriptional or post-transcriptional regulation [24]. Here, we consider the that reduction on RKIP protein numbers happens because of repression of transcription.
Transcription of RKIP can be silenced by methylation of its promoter. Indeed, methylationspecific PCR (MSP) analysis has shown a sufficiently strong correlation between RKIP promoter methylation and low RKIP expression levels in cancer tumors, [24], like esophageal and gastric [34,35]. Histone modifications are also found as epigenetic mechanisms to regulate RKIP levels. Histone deacetylase inhibitors can increase RKIP transcripts [36,37]. Snail and BACH1 transcription factors can downregulate RKIP transcription by histone methyltransferases [38,39] and both are associated with epithelial-mesenchymal transition. Snail is a direct transcriptional repressor of the gene encoding the cell adhesion protein E-cadherin [40], and BACH1 is the basic leucine zipper protein expressed in mammalian tissue which positively regulates motility related genes that promote metastasis in breast cancers [41]. Expression of either Snail and BACH1 genes are self-repressed, and repressed by RKIP proteins. Additionally, BACH1 and RKIP combine into a bistable gene circuit that describe a switch for metastatic phenotype in tumor cellular population, as recently shown [39]. Taken together, those preliminary data indicate that RKIP expression levels may be used as a prognostic marker for survival probabilities and as a target for anti-metastatic treatment.

An effective model for regulation of gene expression
We interpret the gene as a source of gene products randomly switching between ON and OFF states. Synthesis of gene products takes place when the gene is ON at rate k while at OFF state there is no synthesis. The rate of degradation of gene products is denoted by ρ. The gene switches from state ON to OFF, and from OFF to ON, with rates h, and f , respectively.
The aforementioned processes can be represented as a system of effective chemical reactions as given at Eqs. (1)(2)(3)(4). We denote a gene product by P and its regulatory site by R. In this manuscript, we consider the particular case of positive regulation of the gene by a transcription factor denoted by T F . Eqs.
(1-2) indicate, respectively, the gene product synthesis and degradation. The switching from OFF to ON state because of the binding of the activating protein, and the inverse transition caused by its unbinding, are respectively indicated in Eqs. (3)(4). One may also consider the case of a transcription factor being a repressor. Then the effective reactions Eqs. (3)(4) denote the ON to OFF and OFF to ON state transitions, with rates transformed as f → h and h → f . This system of effective reactions is very simple if we consider the complexity of regulation of gene expression in mammals. However, such a simplification enables the construction of exactly solvable quantitative models having a smaller number of parameters which we can use to investigate hypothetical treatment strategies before doing experiments. One example is the use of reduced dose multi-drug treatment targeting a functional network. A given drug has a specific half-life and the multiple drugs act on a system will combine as multi-timescale processes. Combining the dosage and application agenda of those multiple drugs to ensure both effectiveness and non-toxicity may be a hard task and a quantitative model may provide invaluable insights on the therapeutic design. Here we use the stochastic binary model for regulation of gene expression to investigate how the combination of two drugs modulates the dynamics of mRNA production.

An effective model for investigating regulation of gene expression dynamics after treatment
A biological interpretation of the effective model presented at Eqs. (1-4) is given. The rate k is proportional to the inverse of the time interval between two consecutive bindings of RNAPolII to the promoter site. It implies on assuming that transcription starts after a negligible interval following RNAPolII binding. Alternatively, it might be interpreted as the inverse of the average interval between initiation of two subsequent transcription processes. The first case implies on assuming the availability of large amounts of RNAPolII and on interpreting k as a consequence of the affinity between the promoter and RNAPolII. The second could be interpreted as the efficiency of the promoter on initiating transcription. Those two interpretations are not exclusive and, for simplicity, we refer to an increase on k as an increase on the efficiency of the promoter.
The rates of switching h and f are the inverse of the average time of availability and unavailability of the promoter for the binding of the RNAPolII, respectively. The value of those rates will be determined by the binding of transcription factors to the regulatory regions of the gene. Despite the amounts of transcription factors change with time, it is fair to assume that they will remain constant during some sufficiently short intervals. Hence, during those short time intervals, we may employ the effective model of Eqs. (1)(2)(3)(4). The values of the switching constants will reflect the balance resulting from the binding of activators, repressors, quenchers, pioneer factors and other regulatory elements interacting with the regulatory regions of the gene.  Figure 1. Drugs aiming kinetic rates of RKIP transcription. RNAPolII binds to the promoter region (when it is ON) of RKIP gene to synthesize mRNA (gene products). The switching between ON and OFF states is dependent of the regulatory components (proteins) surrounding the gene. Drug 1 targets an activator protein and aims to increase the time of exposition of the promoter for binding of RNAPolII. Drug 2 affects only the promoter by increasing its efficiency or the affinity of the promoter to the RNAPolII. A treatment consists of administering of a single or combinations drugs following a dose agenda.
Our model for regulation of gene expression suggests the design of multiple treatment strategies aiming to increase expression of RKIP gene. The model gives four kinetic rates which we may target by drugs. In the specific case of RKIP gene, one would attempt to increase k and f or to decrease ρ and h. Those rates could be affected individually or collectively, in a coordinated manner, in case of a multi-drug therapy. In the scheme shown in Fig. 1 we consider the particular cases in which the treatment aims to: A) increase f ; B) increase k; C) increase f and k concomitantly. We assume that a given treatment targets a rate exclusively, that is, the non-targeted rates will remain constant during treatment. Two mechanisms can be used for increasing RKIP expression levels: i. Promoter demethylation Downregulation of RKIP has been associated with promoter methylation in many cancer types [34,35,[42][43][44]. We suppose the use of a demethylation agent, namely, 5-Azacytidine (5-AzaC), to revert that (see scheme given in Fig. 2A in [18]). Consequently, expression levels of RKIP would be increased as previously tested in triple-negative breast cancer (TNBC) cell line SUM159 and esophageal cell lines TE-1 and TE-13 [37,42].
ii. Transcription factors regulation The tumor environment is an inexhaustible source of signals which induce a change on the amounts of transcription factors regulating gene expression and, consequently, the phenotypes of tumor cells. For example, nitric oxide (NO) [45] affects tumor growth by downregulating functional quantities of NF-κB and SNAIL which, in reduced quantities, are associated with an increase on RKIP expression. Thus, one may propose a therapeutic strategy using NO donors, such as (Z)-1-[2-(2-aminoethyl)-N-(2ammonio-ethyl) amino] diazen-1-ium-1, 2-diolate (DETANONOate) [40,46,47], to keep a constant inhibition of the NF-κB/SNAIL loop. As a consequence, we expect an increase on the levels of transcripts of RKIP (see scheme in Fig. 2B in [18]).
Finding drugs specifically targeting promoter demethylation of RKIP gene or transcription factors regulating RKIP expression levels is a challenge beyond the scope of this manuscript. Therefore, as a strategy to introduce our methodology, to set a clinically relevant time-scale, and to show the difficulties of combining multiple drugs having multiple targets and time scales we consider non-specific drugs. In that case, we are only considering the effect of the drug on the gene that we describe. The cellular level effects will not be considered here as those would require a more complex approach in which the dynamics of expression levels of multiple potentially interacting genes would need to be considered.

An approach for investigating treatment effects on RKIP expression dynamics
Our treatment strategies aim to increase f and k. We assume that the drug effectiveness on those quantities decays exponentially, and the kinetic rates will return to their pretreatment values. Hence, once treatment starts, the rates f and k become time-dependent with f 0 and k 0 being the OFF to ON and the synthesis rates, respectively, before treatment. We denote the effect of the dose on the value of a kinetic constant by ξ, where 0 < ξ ≤ 1. We assume that at the maximal tolerated dose the effect of the drug is given by ξ = 1 and that that dose raises the targeted rates to f 1 and k 1 . That does not imply on the assumption of maximal efficiency of the drug effectiveness. Indeed, we consider that ξ is the net effect of a given dose on its targeted rate. Hence, a given dose smaller than the maximal tolerated one will instantaneously affect its target rate towards ξ a f 1 and ξ b k 1 , where ξ a and ξ b are non-linear functions of the dose. The time-dependent OFF to ON switching rate, f (t), is and k(t) where j = 1, . . . , J − 1 denotes the j-th drug application, J is the amount of drug doses, τ j is the time of the j-th drug application. At each application of the drug, the steady state condition is that of the untreated system, namely f 0 and k 0 . (λ a , λ b ) denote the rates of exponential decay of the effect of the drugs on the rates ( f , k). As previously mentioned, we are considering that the effect of the drugs on their targeted rates is fast enough to be considered as instantaneous. Therefore, the values of the rates f and k immediately after the arrival of the j-th dose, respectively denoted by f s (τ j ) and k s (τ j ), will be considered as the initial conditions during the interval τ j ≤ t < τ j+1 . Note that the j-th dose adds up to the drug amounts that are remainders from previous applications. Hence, the initial condition at the time τ j is written as: where the dose may be calibrated to generate a differential effect at each time instant to prevent toxic accumulation of drug quantities.

An approximate description of the stochastic binary gene expression dynamics with time-dependent kinetic rates
The randomness of intracellular phenomena suggests a description of the treatment effects on expression of RKIP gene to be built in terms of an stochastic process. Fig 1  suggests the existence of two random variables that determine the state of the system: the gene state, being ON or OFF, and the number of gene products, denoted by n. The description of the dynamics of the state of the system is given in terms of a probability distribution Π(α n (t), β n (t)), where α n (t), or β n (t), denote the probability of finding n proteins at time t when the gene is ON, or OFF, respectively. A master equation governing the probability distribution for an externally regulated gene can be written as: where h and ρ are constants and Eqs. (5) and (6) give f (t) and k(t). The existence of time-dependent coefficients turns hard for one to solve Eqs. (10) and (11) analytically, and some numerical method needs to be employed. The interpretation of the master equation is built in terms of the coefficients k(t), ρ, f (t), and h as presented on the description of the Eqs. (1)(2)(3)(4) and the cartoon shown on Fig. 1.
Here we propose to one to approximate the dynamics of the kinetic rates f (t) and k(t) of Eqs. (5) and (6) as piecewise functions assuming constant values during sufficiently short time intervals ensuring that the difference between the exact and approximated values lies within a given error size. Then, we may consider the model for constant kinetic rates during an interval which is exactly solvable. In that case, the initial condition of an interval is the final condition on its previous neighbor.

An exactly solvable model for benchmarking cancer treatment aiming to modulate gene expression levels
The master equation (10) and (11) with constant coefficients has been already proposed, and it is fully solvable at the stationary [16] and time-dependent regimes [17]. The existence of exact solutions enables one to calculate the time-dependent functions governing the first and the second moment of the number of gene products [18]. Here we write the explicit expressions governing the dynamics of the average number of gene products, n (t) and the standard deviation, σ(t) = n 2 (t) − n 2 (t). We use Eqs. (1)(2)(3)(4) and define the following constants: which are, respectively, the steady state expected number of gene products in case of a gene being fully ON, (we call it maximal mRNA number N); the steady state probability for the gene to be ON (A s ); the ratio of the gene switching rate between ON and OFF states to the degradation rate of the gene products (we call it switching speed ). The average number of mRNAs and the standard deviation at the steady state regime are, respectively, denoted by n s and σ s . We write them as functions of the parameters of Eq. (12): The dynamics of the average and standard deviation is denoted by n (t) and σ(t), respectively, and we write: The coefficients of the exponentials are integration constants given on the Appendix A.2. These solutions are obtained from a system of ordinary differential equation coupling the moments A, n α , and n 2 which exact forms are given on Appendix A.3. Eqs. (15) and (16) enable us to compute the evolution of the average number of products from RKIP gene and its standard deviation using a piecewise representation of the time-dependent rates f (t) and k(t) as we show on next section.
The decaying rate to steady state of both n (t) (Eq. 15) and σ(t) (Eq. 16) can be established in terms of . For 1, the terms e − ρt , e −( +1)ρt become null much faster than the term e −ρt , which hence determines the approaching to steady state. Alternatively, for 1, the term e −ρt will govern lifetime of the dynamical regime, that is, the regime during which n (t) and σ(t) are varying with time.
Additionally, one may notice that also reflects the ratio of the gene switching frequence to the degradation rate, the characteristic times of the two processes being coupled. Eq. (14) indicates that σ 2 s → n s when N > 1. That coincides with the decaying to steady state being determined by the gene product degradation rate, as it happens on a Poisson process. On the other hand, when the gene switching decaying rate to steady state is prevalent, for 1, σ 2 s will have larger values. Then, the gene switching will have a stronger effect on the fluctuations of the number of gene products.

Parameter values and conditions for treatment simulations
We describe simulation conditions for the three treatment scenarios. We denote the degradation rate of DETANONOate (acting on f ) and 5-AzaC (acting on k) by λ a and λ b , which values are λ a = 0.05 h −1 [48] and λ b = 0.25 h −1 (DB00928 entry at DrugBank [49]). The response generated on the kinetic rates by treatment at maximal tolerated dose is indicated by ξ a = ξ b = 1. For a single dose, we assume that the drugs change the values of their targeted rates instantaneously after application, and for maximum tolerated dose, f → f 1 and k → k 1 , as indicated by Eqs. (7) and (8), respectively. Note that the maximum tolerated dose is drug specific and its effect on each kinetic rate can be denoted by the value of ξ.
The gene product of RKIP is assumed to be the mRNAs, which degradation rate is ρ = 0.17 h −1 [50]. The pre-treatment condition, occurring for 0 ≤ t < τ 1 , is not be shown because we assume it as a stationary state so that we set τ 1 = 0. The pre-treatment condition is characterized by a low average copy number of RKIP mRNA's (here chosen as n 0 = 10). We arbitrarily assume that in pre-treatment condition, the levels of mRNA's are 8× below what would be expected to be found in a non-metastatic cell. Hence, the treatment is assumed to be successful if it drives the probability of finding less than n T ≈ 80 mRNAs to a negligible value. That is an important requirement to minimize heterogeneity on treatment response.
The randomness of intracellular processes causes a variability on treatment response even under the hypothetical condition in which all individuals of a population of genetically identical cells absorb the same drug dosage. The heterogeneity of the response can be quantified by the standard deviation of the number of gene products. Hence, we first set n 1 = 100 as the expression level aimed by treatment. This value is chosen assuming that a gene expressing an average of 100 mRNAs in a Poisson regime has the threshold value 80 = n 1 − 2σ 1 . Here, σ 1 = n 1 is the standard deviation of a Poisson distribution having n 1 as its average.
The dynamics of f and k are described by Eqs. (5) and (6). We approximate them as piecewise functions and compute the error using integrals of the exponential decays along each subinterval. Then, we set the length of each subinterval by fixing its error. The piecewise approximation enables the use of n (t) (Eq. 15) and σ(t) (Eq. 16) to describe the expression of RKIP.
A single dose will not be sufficient for keeping n 1 ≈ 100 because of the exponential decay of the drug effect on the kinetic constants (as it will be shown on graphs A-E of Figs. [3][4][5]. Hence, multiple doses are necessary and we determine the intervals between applications of DETANONOate and 5-AzaC as, respectively, 10 h and 4 h. These numbers were chosen to ensure that n (t) ≈ 100 during a sufficiently long time interval and their choice is based on the degradation rates of each drug.
The treatment changes the dynamical properties of the gene switching. When we consider DETANONOate, the pre-treatment probability of finding the gene ON is Then, the aimed steady state probability for the gene to be ON is The preand instantaneously post -treatment value of the gene switching frequence is, respectively, When we consider 5-AzaC the probability of finding the gene ON, given by A 0 = f 0 f 0 +h , and gene switching frequence 0 = f 0 +h ρ , will both remain constant while the value of k 0 → k s (τ j ) instantaneously after drug application at τ j -see

Results
We simulate three treatment scenarios: 1. denotes the effect of maximum tolerated dose of DETANONOate; 2. denotes the effect of maximum tolerated dose of 5-AzaC; 3. denotes the effect of fractionary doses of the two drugs together.
The treatment response is quantified in terms of the time for the average value to reach the threshold and the standard deviation. These two quantities are strongly dependent of the values of 0 , as we will show on next subsections. The values of ρ, λ a , and λ b are assumed to remain constant during treatment. All rates (k, f , h) are given in h −1 .
The values of the parameters are selected for simulating treatment conditions which probability distributions governing expression of pre-treated cells indicate qualitatively distinguishable steady state regimens as recently classified [19]. Fig. 2 shows pre-treatment (and aimed post-treatment) steady state probability distributions governing the numbers of RKIP mRNAs in red (and green). The values of the parameters in each row are set to ensure that graphs A indicate bimodal distributions, graphs B are distributions close to the limit of the bimodal regime, graphs C indicate the regime in which the probabilities are table-shaped for A 0 = A 1 = 0.5 (as indicated in graph C2), graphs D indicate the quasi-Poisson distributions, and graphs E denote the bursting limit. The curves of the steady state probability distributions of finding n mRNAs, within the Fig. 2, are computed by Eq. (A1) and it is denoted byφ n in the Appendix A.1.
The parameters of the pre-and post-treatment probability distributions are fixed such that the steady state average number of mRNA's will be ∼ 10 and ∼ 100, respectively. The bimodal distributions indicate that the mRNAs synthesized while the gene is ON will degrade fastly after the switching to the OFF state. For the table-shaped and quasi- Poissonian limits, the switching between ON-OFF states are, respectively, slow and fast enough to ensure that the mRNA degradation is compensated by its synthesis to generate the specific distributions. The burst limit is characterized by very short ON states during which the synthesis is very efficient while the OFF states duration is proportionally very long. On next subsections we present the results of simulations of the treatment response dynamics considering the five pre-treatment conditions shown in Fig. 2  The first and second rows of graphs of Fig. 3 show the simulation of the dynamics of 1 expression after one and five drug doses, respectively, under 5 pre-treatment conditions. this case the average number does not cross the threshold and the lower line n (t) − σ(t) 12 does not appear. 13 The response to multi-doses treatment is shown in graphs F1 to J1. The interval 14 between doses was chosen to enable a sigmoidal-like response characterizing two distinct which establishes the response to treatment having the minimal heterogeneity. Graph J1 22 shows the response of a bursty gene which n (t) crosses the threshold after ∼ 12 h and 23 reaches a maximal value of ∼ 170. However, this regime causes the noisiest response 24 as indicated by n (t) − σ(t) not crossing the threshold. Note also that there are some 25 "bumps" because the effect of the single dose is close to the maximum by the time of the 26 next drug application. It is possible to prevent the bumps but we decided to show them to 27 demonstrate that the time interval between drug dosages also requires one to consider the 28 dynamics of the targeted gene.  The first and second rows of graphs of Fig. 4 show the simulation of the dynamics of 38 expression after one and ten drug doses, respectively, under 5 pre-treatment conditions. 39 The treatment do not affect either or ρ, the decaying rates of the system are not n (t) − σ(t) remains negative while n (t) + σ(t) exceeds the threshold more than 2×. 46 The response to multi-doses treatment is shown in graphs F2 to J2. The interval 47 between doses was chosen to enable a sigmoidal-like response with the average number 48 reaching at least 100. All graphs show that n (t) cross the threshold after ∼ 10 h. That    Table 1.

71
Graph A3 shows the dynamics of treatment response when the initial condition is set 72 for a slow switching gene. The average number reaches a maximum of ∼ 25 at t ∼ 35 h.

73
The maximum of n (t) + σ(t) is ∼ 70, a highly heterogeneous response if we consider that 74 the maximum standard deviation is equal to the maximum average, that is, it is a super- pre-treatment conditions, the average number has not crossed the threshold.

83
The response to multi-doses treatment is shown in graphs F3 to J3. Table 1 shows the 84 sequences of applications of a given drug in fractions of maximal tolerated dose, ξ a and ξ b .

85
The values were chosen aiming that the total amount of drugs is reduced in comparison 86 with the single drug treatment. We also attempt to ensure a sigmoidal-like response such 87 that the average number will reach at least ∼ 100. Graph F3 shows the dynamics of  Table 1). The agenda enables a reduction of 20% in comparison with application 90 of full doses (ξ a = ξ b = 1). Because the pre-treatment gene is in a slow switching regime 91 that leads to the slowest response as n (t) crosses the threshold after ∼ 40 h. It also has 92 Table 1. The fractions ξ a and ξ b of maximal tolerated dose for treatment agendas of pre-treatment conditions in Fig. 5 Graph Sequence of number of doses × fraction for Cumulative reduction in a high noise as it can be noticed by the values of n (t) ± σ(t). Graphs G to J3 show the 93 simulations of the response for 6 (or 15) doses of drugs a (or b) with ξ a and ξ b ranging from 94 0.8 to 0.5 (see the Table 1). The cumulative reduction in full doses range from 29% to 35% 95 (or from 23% to 35%) in ξ a (or ξ b ). Graphs G3 to I3 show that n (t) crosses the threshold 96 after ∼ 15 h. The curves n (t) ± σ(t) become closer as the value of increases which 97 indicates the direction to design treatment strategies with reduced response heterogeneity.

98
Graph J3 shows the bursty pre-treatment condition. Here, n (t) crosses the threshold after 99 ∼ 22 h and n (t) + σ(t) reaches ∼ 200, which indicates the noisiest response. set based on mRNA-microRNA binding duration time [51] and h in BACH1 half-life [52].

115
The time-dependence of h(t) and ρ(t) is described following the same framework used for 116 k(t) and h(t) -see Eqs. (5)(6)(7)(8). To obtain the dynamics of the average number of mRNAs 117 and of the standard deviation, we extend the previous formulation of the master equation 118 of Eqs. (10) and (11) to all kinetic rates.

119
The standard deviation of the post-treatment probability distribution is set to be which has a local minimum at N 1 → 0 and a local maximum at where N is the maximal mRNA number, A is the steady state probability for the gene 120 to be ON, is the gene switching speed (Eq. 12) and the subscript 1 indicates aimed 121 post-treatment parameter values. A strategy to reduce σ 1 is to keep N 1 1 + 1 (or 122 k 1 f 1 + h 1 + ρ 1 ), which ensures that σ 1 approaches the local maximum as A 1 → 1. The 123 reduction of σ 1 values depends on the mean number n 1 desired, so that n 1 ≈ A 1 (1 + 1 ).

154
When the binary model for regulation of gene expression is considered, it is useful to 155 characterize the pre-treatment gene expression regime in terms of A 0 , 0 , and N 0 (see Eq. 156 12). Those rates will affect the pre-existing (and post-treatment) noise on the number of 157 mRNAs, see Eqs. 13 and 14 (and Figs. 3,4,5 and 6). Furthermore, the time of response to 158 treatment also depends on 0 , because the decaying to steady state regime of the average 159 number of mRNAs and its variance have both ρ and ρ as their smallest rates (see Eq. 15).

162
All trajectories of the average values shown in graphs of Fig. 4  despite treatment increasing f will cause an increase of the relative switching speed.

193
In the high noise responses the trajectories n (t) − σ(t), do not cross the threshold. It  The synthesis of RKIP protein may have a non-linear dependence of the mRNA 207 numbers, because translation can be regulated by mechanisms such as mRNA stability, 208 translational control and proteosomal degradation [50]. Hence, a theoretical approach 209 considering translation might required the addition of two or more characteristic timescales 210 depending on protein synthesis and degradation [53]. We emphasize that the amounts 211 of RKIP mRNAs in our model are the net effect of the enhanced treatment, since its 212 expression results from the interaction with a complex gene network. Our approach, 213 however, provides the building blocks of gene networks [15] that can be used to understand 214 the functioning of larger modules [10][11][12] and, hence, to further enhance treatment designs 215 of cancer at metastatic stage.

216
To further improve our approach one needs to perform a data-based validation of 217 our model. One possible experiment would involve the decrease of RKIP levels, using a 218 siRNA approach to simulate different levels of activation of the RKIP gene. This would 219 be a more specific and precise way to interfere in the RKIP pathway than using 5-AzaC. transcripts and promoter activity [28]. This regime also coincides with condition of highly 232 heterogeneous and slow response to treatment. Hence, we may ask whether it is relevant 233 for a master regulatory gene to transmit information about its own promoter state and in 234 which cellular context, namely, in a cancerous or healthy one. The answer to this question 235 is important because it helps us to understand the role of noise and the conditions under 236 which its reduction is desirable [27] as it happens when we have a negatively self-regulating 237 gene. iii. both rates of a gene using both drugs together. We show that treatment response 249 speed and heterogeneity depends on the pre-treatment state of the gene. Then, we present 250 an enhanced treatment design that ensures reduced heterogeneity and time of response.

251
Besides useful to inspect treatment designs, the response to treatment may be used on 252 inference of the kinetic constants of a given gene in a synthetic system.  n α (t) = n α s + Re −(h+ f )t + Se −ρ t + Te −( f +h+ρ)t , (A11) n 2 (t) = n 2 s + Ue −(h+ f )t + (1 + 2 n s )Ve −ρ t + We −( f +h+ρ)t + Xe −2 ρ t , (A12) The following table shows the treatment parameters used to obtain the graphs of