Next Article in Journal
Neoadjuvant Chemotherapy of Patients with Early Breast Cancer Is Associated with Increased Detection of Disseminated Tumor Cells in the Bone Marrow
Next Article in Special Issue
Harnessing RKIP to Combat Heart Disease and Cancer
Previous Article in Journal
Clinical Value of Surveillance 18F-fluorodeoxyglucose PET/CT for Detecting Unsuspected Recurrence or Second Primary Cancer in Non-Small Cell Lung Cancer after Curative Therapy
Previous Article in Special Issue
Insights of RKIP-Derived Suppression of Prostate Cancer
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Stochastic Binary Model for the Regulation of Gene Expression to Investigate Responses to Gene Therapy

by
Guilherme Giovanini
1,
Luciana R. C. Barros
2,
Leonardo R. Gama
2,
Tharcisio C. Tortelli, Jr.
2 and
Alexandre F. Ramos
1,2,*
1
Escola de Artes, Ciências e Humanidades, Universidade de São Paulo, Av. Arlindo Béttio, 1000, São Paulo 03828-000, SP, Brazil
2
Centro de Investigação Translacional em Oncologia, Departamento de Radiologia e Oncologia, Faculdade de Medicina da Universidade de São Paulo, Instituto do Câncer do Estado de São Paulo, Av. Dr. Arnaldo, 251, São Paulo 01246-000, SP, Brazil
*
Author to whom correspondence should be addressed.
Cancers 2022, 14(3), 633; https://doi.org/10.3390/cancers14030633
Submission received: 7 September 2021 / Revised: 8 November 2021 / Accepted: 13 November 2021 / Published: 27 January 2022
(This article belongs to the Special Issue RKIP: A Pivotal Gene Product in the Pathogenesis of Cancer)

Abstract

:

Simple Summary

Gene editing technologies reached a turning point toward epigenetic modulation for cancer treatment. Gene networks are complex systems composed of multiple non-trivially coupled elements capable of reliably processing dynamical information from the environment despite unavoidable randomness. However, this functionality is lost when the cells are in a diseased state. Hence, gene-editing-based therapeutic design can be viewed as a gene network dynamics modulation toward a healthy state. Enhancement of this control relies on mathematical models capable of effectively describing the regulation of stochastic gene expression. We use a two-state stochastic model for gene expression to investigate treatment response with a switching target gene. We show the necessity of modulating multiple gene-expression-related processes to reach a heterogeneity-reduced specific response using epigenetic-targeting cancer treatment designs. Our approach can be used as an additional tool for developing epigenetic-targeting treatments.

Abstract

In this manuscript, we use an exactly solvable stochastic binary model for the regulation of gene expression to analyze the dynamics of response to a treatment aiming to modulate the number of transcripts of a master regulatory switching gene. The challenge is to combine multiple processes with different time scales to control the treatment response by a switching gene in an unavoidable noisy environment. To establish biologically relevant timescales for the parameters of the model, we select the RKIP gene and two non-specific drugs already known for changing RKIP levels in cancer cells. We demonstrate the usefulness of our method simulating three treatment scenarios aiming to reestablish RKIP gene expression dynamics toward a pre-cancerous state: (1) to increase the promoter’s ON state duration; (2) to increase the mRNAs’ synthesis rate; and (3) 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 reaching increased average mRNA levels with diminished heterogeneity while reducing drug dosage by simultaneously targeting multiple kinetic rates that effectively represent the chemical processes underlying the regulation of gene expression. The decrease in heterogeneity of treatment response by a target gene helps to lower the chances of emergence of resistance. Our approach may be useful for inferring kinetic constants related to the expression of antimetastatic genes or oncogenes and for the design of multi-drug therapeutic strategies targeting the processes underpinning the expression of master regulatory genes.

1. Introduction

Recent advances in gene editing technologies brought the promise of a turning point for gene therapy [1] toward more complex therapeutic designs aiming to orchestrate the expression of gene networks for cell phenotype reprogramming [2]. One possibility is to develop cancer treatment strategies to revert metastasis by targeting master regulatory genes [3]. Mathematical models describing the regulation of gene expression can be insightful for engineering of the dynamics of the gene networks governing cellular behavior.
Let us assume the ideal case in which the editing exclusively affects its epigenetic target [2] within tumor cells. The task can be formulated as a control problem to enable the number of transcripts of a master regulating gene to have its average value at a given level and random fluctuations within a sufficiently small range. The control may be performed by external agents, such as a combination of drugs that we would like to keep at a minimally effective dosage because of the eventual toxicity.
Despite our deepened understanding of cancer biology due to the advances in molecular biology techniques, the use of quantitative methods to integrate the plethora of generated data to design treatments targeting metastasis is still in its infancy [4,5]. The genotypic variability and intrinsic randomness of biochemical reactions [6] governing epigenetics underpin the multiplicity of cancer cell phenotypes commonly termed as tumor heterogeneity [7,8,9,10,11,12].
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 [13,14,15,16]. The numerous characteristic timescales of the chemical processes taking place inside the cell add another layer of cumbersomeness [17,18].
Thus, the enhancement of therapeutic strategies requires the analysis and engineering of the dynamics of treatment response in a complex system composed of several interacting components with a multiplicity of characteristic time scales and subjected to randomness. Finding building blocks with those features may provide useful insights on how to modulate the dynamics of such a complex system [19]. For that, the exactly solvable stochastic model for transcription of a binary gene [20,21,22,23] is a good candidate to be used as a prototype for simulating enhanced treatment strategies.

Rationale

In this manuscript, we propose a proof of concept of a quantitative tool to investigate the response to the application of multiple drugs for epigenetic control of master regulatory genes. Those treatment designs involve a large number of chemical compounds with diverse half-lives and interacting affinities with the gene components. A gene may have its transcription guided by a promoter with multiple states, and the duration of those states may be regulated by multiple transcription factors with various affinities to the regulatory sites of the gene.
For the promoter in the ON state, one would still have variation on the affinity between PolII and the promoter. Additionally, one needs to select the typical external agents, here drugs, that will affect the processes involved in the modulation of gene expression accordingly with their pharmacokinetic parameters. Hence, mathematical modeling of treatment targeting specific genes requires the selection of a proper gene model system to enable the determination of the parameter values associated with the regulation of its expression. One needs to couple this model with the pharmacokinetics of the interaction of the drugs with their targets.
This can be formulated as a control problem in which we aim to keep both the average expression level of a target gene and the heterogeneity of its response to treatment within specific ranges while minimizing the amounts of drugs used to prevent toxicity. To devise a quantitative strategy to overcome the challenge above, we started with the simplest possible mathematical model for regulated gene expression under influence of randomness to establish a proof of concept. This stochastic binary model for gene expression has four parameters and two characteristic time scales. Ideally, each of the multiple drugs composing a treatment would affect a specific process (or mathematically speaking, parameter) involved in the expression of the target gene.
This implies on adding at least two parameters per drug: the halflife of the drug, and the strength of its effect in each specific parameter of the gene expression model. The next challenge is to orchestrate this multiplicity of parameters to properly regulate the expression of the target gene. Here, we analyze the drug dose and gene response using a “pedestrian optimization” as application of optimal control theory is a highly elaborate technique with specific requirements and the formulation is beyond the scope of the current study [24]. One advantage is that multiple qualitative features of the stochastic model for gene expression are preserved independently of the numerical values that are used.
Our next challenge is the selection of an appropriate candidate model system. Here, we use the gene encoding Raf kinase inhibitory protein (RKIP) because it plays a major role in regulating the dynamics of multiple components of a cell [25]. Since cancer lethality is mainly caused by metastasis, the choice of RKIP is promising as its concentrations are typically reduced in metastatic cancers [26,27,28]. This negative correlation turns RKIP and RKIP-related gene signatures into useful biomarkers of metastatic risk in cancer patients [29].
The characterization of RKIP as an anti-metastatic gene is reinforced by experiments demonstrating that its overexpression blocks in vivo invasion and metastatic progression [30,31]. Hence, we focus on possible strategies to increase the amounts of RKIP in cancer cells by re-modulating its expression profiles toward the pre-cancer regimens. In some cancers, the mean number of RKIP transcripts is similar to that of cells of a normal tissue, but the variability is significantly increased [32]. Our approach is also potentially useful for those cases because it is based on stochastic processes.
The selection of RKIP as a model system provides us with its degradation rate, i.e., one of the characteristic times of our stochastic model for regulation of gene expression. The additional characteristic time is the gene switching rate, which is usually unavailable because it needs to be inferred. Once the RKIP gene is selected, we may also consider a multi-drug treatment based on 5-AzaC and DETANONOate because they have been tested previously, their mechanisms of action are sufficiently known, and because we could recover their half-life times.
Hence, in the case of a treatment design with two drugs, we will have three out of four half-life times for the system composed by the genes and treatment. This prompts us to investigate the qualitative features of our prototypic model as an additional tool for simulating cancer treatment designs targeting specific genes (Note, however, that parameter value adjustments for specific experimental designs can be performed if one judges that our model may provide sufficiently useful results).
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. Due to the treatment, we will consider the kinetic rates of the model to be time-dependent in response to drug application. This 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 [23,33] and the reliability of information about the promoter state that is transmitted by gene products [34]. Here, we assume that the addition of a drug causes the kinetic rates of the model to become time dependent and that this effect decays exponentially such that it has up to four additional time scales for the problem. Then, we show that treatment targeting a gene with a fast switching promoter, and expressed 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 expressed in a burst fashion will enable a fast response with 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 [4,5].
The remainder of this manuscript is organized as follows. The assumptions underlying the qualitative and quantitative models that we propose are presented in Section 2. The obtained results are shown in Section 3, and we discuss them in Section 4. Our concluding remarks are presented in Section 5.

2. Methods and Models

Gene expression is the key mechanism in cell function by means of an extensive molecular machinery that transforms the genetic information into molecular functions. Gene expression can be described as a two-stage process—namely, the gene transcription that generates mRNA, and the mRNA translation, which has proteins as products. Here, we focus on gene transcription and provide a simplified picture of how it occurs 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 an enhancer. The enhancers are regions of 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.

2.1. A Brief Description of the Molecular Role of RKIP

RKIP protein is a regulator of kinases that directly binds to Raf kinase [32,35,36,37]. RKIP is involved in regulation of signaling pathways, such as the Raf–Mek–Erk cascade and the NF- κ B-related pathways [38,39], both participating in the regulation of anti-apoptosis processes. Additionally, those pathways modulate cell proliferation with the Raf–Mek–Erk cascade participating in differentiation and NF- κ B-related pathways acting in inflammation. In the Raf–Mek–Erk cascade, RKIP inhibits the downstream signaling pathway by the direct binding of the dephosphorylated RKIP protein to Raf-1.
This molecular complex prevents both Raf-1 phosphorylation and Raf-1/Mek association, which, in turn, causes the interruption of Erk signaling. This inhibition can be reverted by action of Protein Kinase C (PKC) [40,41,42], a post-translational regulator that phosphorylates RKIP. The latter becomes dissociated from Raf-1 because of its structural change, and the Erk pathway becomes activated. RKIP negatively regulates the NF- κ B signaling pathway indirectly [32,38] when it interacts with kinase complex IKK reducing their activity.
This causes a reduction in phosphorylation and degradation of the inhibitory proteins I κ B, which, in turn, inactivate NF- κ B. In addition, RKIP modulates other fundamental cell signaling pathways involving heterotrimeric G-proteins, keap1/nrf2, STAT3, and GSK [43]. Due to 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 [44]. A plethora of metastatic solid tumors have RKIP downregulated or lost, and the experimental data suggests that this happens because of transcriptional or post-transcriptional regulation [29]. Here, we limit our attention to cases in which a reduction in RKIP protein numbers happens because of repression of transcription.
Transcription of RKIP can be silenced by methylation of its promoter. Indeed, methylation-specific PCR (MSP) analysis has shown a sufficiently strong correlation between RKIP-promoter methylation and low RKIP expression levels in cancer tumors, [29], including esophageal and gastric [45,46]. Histone modifications are also found as epigenetic mechanisms to regulate RKIP levels.
Histone deacetylase inhibitors can increase RKIP transcripts [47,48]. Snail and BACH1 transcription factors can downregulate RKIP transcription by histone methyltransferases [49,50], and both are associated with the epithelial–mesenchymal transition. Snail is a direct transcriptional repressor of the gene encoding the cell adhesion protein E-cadherin [51], and BACH1 is the basic leucine zipper protein expressed in mammalian tissue, which positively regulates motility related genes that promote metastasis in breast cancers [52]. The 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 describes a switch for metastatic phenotype in tumor cellular population, as recently shown [50].
RKIP mRNAs may also be post-transcriptionally regulated negatively because of the interaction with specific miRNAs. Indeed, in several cancers, the RKIP gene is silenced by the action of miRNAs (-224, -27a, -23a, and -543) targeting RKIP mRNAs [53,54]. This suggests a therapeutic alternative based on targeting miRNAs that downregulate RKIP expression, such as lncRNA XIST, which stabilizes RKIP expression by suppressing miR-23a [55].
The aforementioned data suggests an association between RKIP expression levels and tumor cell phenotypes. This leads to the possibility of using RKIP as a prognostic marker for survival probabilities [56,57]. Additionally, as RKIP levels are mainly reduced in several cancers [32], understanding its regulatory mechanisms [45,46,47,48,49,50,51,52] may enable the design of anti-metastatic treatments.

2.2. An Effective Model for the 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 Equations (1)–(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 .
RT F k RT F + P ,
P ρ ,
R + T F f RT F ,
RT F h R + T F .
Equations (1) and (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 Equations (3) and (4). One may also consider the case of a transcription factor as a repressor. Then, the effective reactions Equations (3) and (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 the regulation of gene expression in mammals. However, such a simplification enables the construction of exactly solvable quantitative models with a smaller number of parameters that we can use to investigate hypothetical treatment strategies before performing 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 that act on a system will combine as multi-timescale processes.
Combining the dosage and application agendas 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 the regulation of gene expression to investigate how the combination of two drugs modulates the dynamics of mRNA production.

2.3. An Effective Model for Investigating Regulation of Gene Expression Dynamics after Treatment

A biological interpretation of the effective model presented at Equations (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. This implies assuming that transcription starts after a negligible interval following RNAPolII binding. Alternatively, it might be interpreted as the inverse of the average interval between the initiation of two subsequent transcription processes.
The first case implies assuming the availability of large amounts of RNAPolII and 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 in k as an increase in 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 amount of transcription factors that change with time, it is fair to assume that they will remain constant during sufficiently short intervals. Hence, during those short time intervals, we may employ the effective model of Equations (1)–(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.
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 that we may target by drugs. In the specific case of the RKIP gene, one can attempt to increase k and f or to decrease ρ and h. Those rates can be affected individually or collectively, in a coordinated manner, in the case of a multi-drug therapy. In the scheme shown in Figure 1, we consider the particular cases in which the treatment aims to: (1) increase f; (2) increase k; and (3) 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:
(1) Promoter demethylation. The downregulation of RKIP has been associated with promoter methylation in many cancer types [45,46,58,59,60]. We suppose the use of a demethylation agent, namely, 5-Azacytidine (5-AzaC), to revert that (see scheme given in Figure 2A in [22]). Consequently, the expression levels of RKIP would be increased as previously tested in the triple-negative breast cancer (TNBC) cell line SUM159 and esophageal cell lines TE-1 and TE-13 [48,58].
(2) Transcription factors regulation. The tumor environment is an inexhaustible source of signals that induce a change in the amounts of transcription factors regulating gene expression and, consequently, the phenotypes of tumor cells. One possible treatment is to use nitric oxide (NO) or NO donors acting as direct or chemosensitizer anti-cancer agents [61]. For example, NO [62] affects tumor growth by downregulating the functional quantities of NF- κ B and SNAIL, which, in reduced quantities, are associated with an increase in RKIP expression.
Thus, one may propose a therapeutic strategy using NO donors, such as (Z)-1-[2-(2-aminoethyl)-N-(2-ammonio-ethyl) amino] diazen-1-ium-1, 2-diolate (DETANONOate) [32,51,63,64], to maintain a constant inhibition of the NF- κ B/SNAIL loop. Indeed, RKIP was increased in treatments with DETANONOate designed for inhibiting the epithelial–mesenchymal transition (EMT) and invasion in metastatic human prostate carcinoma cell lines, which was corroborated in mice bearing tumor xenografts [51].
NO donors were used in combination with photosensitizers to increase the efficacy of the photodynamic therapy inhibiting proliferation of murine melanoma cells [65]. Note, however, that NO plays a dual role [66] since, in low-doses, it promotes carcinogenesis [67,68] Indeed, low levels of photodynamic therapy induce low levels of NO, which contributes to an anti-apoptotic response by NF- κ B/Snail/RKIP loop [69]. In this manuscript, we assume that the concentrations of NO within cancer cells are sufficiently large (micromolar levels [70]) to ensure its role in promoting an increase in the amount of RKIP transcripts (see scheme in Figure 2B in [22]).
The choice of 5-AzaC and DETANONOate was motivated because of a sufficiently good characterization of their mechanisms of action on signaling pathways that participate in upregulating RKIP expression. However, alternative treatment designs might be used to promote RKIP expression enabling a synergistic association with cancer therapies, such as Sorafenib associated with Gemcitibine [71] or Erlotinib [72] in lung cancers, and Gemcitibine with Sorafenib in pancreatic cancer [73]. Topoisomerase I inhibitor 9NC [74] and anti-mitotic agents ENMD-1198 and MKC-1 have also demonstrated the upregulating effects of RKIP [75].
However, finding drugs specifically targeting the 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 clinically relevant timescales, and to show the difficulties of combining multiple drugs with multiple targets and timescales, 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 the expression levels of multiple potentially interacting genes would need to be considered.

2.4. 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 that the kinetic rates will return to their pre-treatment 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 fractional effect of a drug dose on the value of a kinetic constant by ξ , where 0 < ξ 1 . We assume that, at the maximum tolerated dose, the effect of the drug is given by ξ = 1 and that this dose raises the targeted rates to f 1 and k 1 .
This does not imply 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 maximum tolerated one will instantaneously change its target rate to ξ a f 1 and ξ b k 1 , where ξ a and ξ b are non-linear functions of the dose that need to be formulated accordingly with experimental data. The time-dependent OFF to ON switching rate, f ( t ) , is
f ( t ) = f 0 , 0 t < τ 1 f 0 + f s ( τ j ) f 0 e λ a ( t τ j ) , τ j t < τ j + 1
and k ( t )
k ( t ) = k 0 , 0 t < τ 1 k 0 + k s ( τ j ) k 0 e λ b ( t τ j ) , τ j t < τ j + 1
where j = 1 , , J 1 denotes the j-th drug application, J is the amount of drug doses, and τ 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:
f s ( τ j ) = f 1 i = 1 j ξ a ( τ i ) e λ a ( τ j τ i ) ,
k s ( τ j ) = k 1 i = 1 j ξ b ( τ i ) e λ b ( τ j τ i ) ,
where the dose may be calibrated to generate a differential effect at each instant to prevent toxic accumulation of drug quantities.

2.5. 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 a stochastic process. Figure 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:
d α n d t = k ( t ) ( α n 1 α n ) + ρ [ ( n + 1 ) α n + 1 n α n ] h α n + f ( t ) β n ,
d β n d t = ρ [ ( n + 1 ) β n + 1 n β n ] + h α n f ( t ) β n ,
where h and ρ are constants and Equations (5) and (6) give f ( t ) and k ( t ) . The existence of time-dependent coefficients is difficult to solve Equations (10) and (11) analytically, and some numerical methods need 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 Equations (1)–(4) and the cartoon shown on Figure 1.
Here, we propose to approximate the dynamics of the kinetic rates f ( t ) and k ( t ) of Equations (5) and (6) as piece-wise 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 (see Equation (A1) in the Appendix A.1). Then, we may consider the model for constant kinetic rates during an interval that is exactly solvable. In that case, the initial condition of an interval is the final condition on its previous neighbor.

2.6. An Exactly Solvable Model for Benchmarking Cancer Treatment Aiming to Modulate Gene Expression Levels

The master Equations (10) and (11) with constant coefficients has been already proposed, and it is fully solvable at the stationary [20] and time-dependent regimes [21]. 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 [22]. 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 Equations (1)–(4) and define the following constants:
N = k ρ ; A s = f f + h ; ϵ = f + h ρ ,
which are, respectively, the steady state expected number of gene products in the case of a gene being fully ON, (we call this the maximal mRNA number N); the steady state probability for the gene to be ON ( A s ); and the ratio of the gene switching rate between ON and OFF states to the degradation rate of the gene products (we call this the 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 Equation (12):
n s = A s N ,
σ s 2 = n s 1 + N 1 A s 1 + ϵ .
The dynamics of the average and standard deviation are denoted by n ( t ) and σ ( t ) , respectively, and we write:
n ( t ) = n s + Y e ϵ ρ t + V e ρ t ,
σ 2 ( t ) = σ s 2 + U 1 e ϵ ρ t + V e ρ t + W 1 e ( 1 + ϵ ) ρ t + X 1 e 2 ρ t Y 2 e 2 ϵ ρ t .
The coefficients of the exponentials are integration constants given on the Appendix A.4. These solutions are obtained from a system of ordinary differential equation coupling the moments A, n α , and n 2 , where the exact forms are given in Appendix A.5.
Equations (15) and (16) enable us to compute the evolution of the average number of products from the RKIP gene and its standard deviation using a piece-wise representation of the time-dependent rates f ( t ) and k ( t ) as we show in the next section.
The decaying rate to steady state of both n ( t ) (Equation (15)) and σ ( t ) (Equation (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 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 frequency to the degradation rate, the characteristic times of the two processes being coupled. Equation (14) indicates that σ s 2 n s when ϵ N > 1 . This 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 , σ s 2 will have larger values. Then, the gene switching will have a stronger effect on the fluctuations of the number of gene products.

2.7. Parameter Values and Conditions for Treatment Simulations

The exactly solvable model for binary stochastic gene expression enables the simulation of the dynamics of the response to treatment measured in terms of both the average expression and random fluctuations around the mean of a hypothetical (abstract) gene. However, such a choice may turn hard the provision of insights for one willing to design a treatment to modulate the expression of a master regulating gene. One goal is the design of an optimal agenda of application of treatment to control both the average and standard deviations of gene products to remain within specific ranges.
The instants of application of each drug, denoted by τ , are strongly related to the characteristic times of the system composed of a gene and its interacting drugs. Among those characteristic times, the half-lives of the drugs, denoted by λ i where i labels a given drug, and degradation rate of the gene products, denoted by ρ , are widely used. However, as illustrated in a previous subsection, the introduction of the promoter switching provides an additional characteristic time, denoted by ϵ ρ , which will have very important effects on treatment response, as shown below. Hence, the choice of the system RKIP, 5-AzaC, and DETANONOate helps us to set a few biologically reasonable parameter values to perform our analysis. Adapting our framework for different systems, however, is a fully feasible task.
Here, the rates related to DETANONOate and 5-AzaC are, respectively, denoted by an index a and b. Hence, the degradation rates of DETANONOate (acting on f) and 5-AzaC (acting on k) are denoted by λ a and λ b , respectively. Their values are λ a = 0.05 h 1 [76] and λ b = 0.25 h 1 (DB00928 entry at DrugBank [77]).
A separate challenge for designing a therapy targeting a specific gene is the formulation of an effective model for the effect of treatment on a specific gene parameter. As this is beyond the scope of the current study, we defined a quantity denoted by ξ , which will indicate the effect of a given treatment dose on its target parameter. For the maximum tolerated dose, the effect will be considered as maximal and indicated by ξ = 1 . When we have 0 < ξ < 1 , we are indicating that the dosage is smaller than maximal. Note, however, that we do not have a mathematical model relating ξ to the drug dosage, and its construction requires specific experiments.
We denote the response generated on the kinetic rates by treatment at maximum tolerated dose 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 Equations (7) and (8), respectively.
The specific values of dosage of both drugs considered here have been set in previous studies. The dose values and effects on RKIP expression levels may vary by up to one order of magnitude. For example, experimental analysis of human prostate metastatic cells (DU145 and PC-3) treated with 1000 μ M DETANONOate showed upregulation of RKIP mRNA for 4 and 12 h post-treatment [51]. Triple-negative breast cancer cells (SUM159) also increase in RKIP mRNA 1.4 fold when treated with 500–2000 μ M 5-AzaC for 72 h [48], and human esophageal cancer cell (TE-1 and TE-13 cells) treated with 2 μ g/mL 5-AzaC showed a 1.5 –10-fold increase in RKIP expression [58].
The treatment induces a change in the kinetic parameters of the model. As a consequence, all quantities depending on those parameters will be changed. Hence, f 0 f 1 and k 0 k 1 causes a change of ϵ 0 ϵ 1 , and a change on the steady state values of the statistical quantities n s , 0 n s , 1 , and σ s , 0 σ s , 1 accordingly with Equations (12)–(14). The expected steady state probability distributions governing the number of RNA transcripts instantaneously after treatment application as shown in Figure 2. However, because the drug effects decay exponentially with a rate determined by their half-life, the system does not reach those new steady state values and tends to return to the pre-treatment conditions instead.
Here, we assume that the gene product of RKIP are mRNAs, whose degradation rate is ρ = 0.17 h 1 [78]. The pre-treatment condition, occurring for 0 t < τ 1 , is not shown because we assume it as a stationary state, and thus 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 conditions, the levels of mRNAs are eight-times 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. This is an important requirement to minimize heterogeneity in a treatment response.
The randomness of intracellular processes causes a variability in the treatment response even under the hypothetical conditions 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 with n 1 as its average.
The dynamics of f and k are described by Equations (5) and (6). We approximate them as piece-wise 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 piece-wise approximation enables the use of n ( t ) (Equation (15)) and σ ( t ) (Equation (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 will be shown on graphs A–E of Figure 3, Figure 4 and Figure 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 A 0 = f 0 f 0 + h . The drug dose delivered at τ j causes f 0 f s ( τ j ) instantaneously as defined at Equation (7). Then, the aimed steady state probability for the gene to be ON is A 1 = f s ( τ j ) f s ( τ j ) + h . The pre (and instantaneously post) treatment values of the gene switching frequency are, respectively, ϵ 0 = f 0 + h ρ and ϵ 1 = f s ( τ j ) + h ρ . When we consider 5-AzaC the probability of finding the gene ON, given by A 0 = f 0 f 0 + h , and gene switching frequency ϵ 0 = f 0 + h ρ , will both remain constant, while the value of k 0 k s ( τ j ) instantaneously after drug application at τ j —see Equation (8)—such that N 0 = k ρ N 1 = k s ( τ j ) ρ .
As we are using an exactly solvable stochastic model, we can map its qualitative features in terms of the relations between its kinetic parameters. For example, the value of ϵ in Equation (12) is a key parameter determining the shape of the steady state probability distributions shown in Figure 2. Hence, the interpretation of our results will remain useful in the analysis of specific experimental designs despite the need for specific values for the parameters of the model.
Hence, our choice for both pre- and post-treatment average values of mRNAs are arbitrary and selected for the clearness of presentation of our results and predictions. A four- to eight-fold increase in the median of mRNA numbers has been reported in PCPG (pheochromocytoma and paraganglioma), CHOL (cholangiocarcinoma), and SARC (sarcoma) as shown in [32], Figure 1. The RKIP expression levels in LIHC (liver hepatocellular carcinoma) [79], GBMLGG (glioma) [80] and STES (stomach and esophageal carcinoma) [81] were increased by up to two times.
Other comparisons of RKIP expression levels in breast cancer (BRCA and TNBC) [82,83], skin cancer (SKCM) [84,85], and colorectal cancers (COADREAD) [57,86,87] indicated a negligible change on the median. In all aforementioned cancers, the variability in the numbers of transcripts of RKIP was larger in tumor cells than in normal cells. One advantage of the use of a stochastic model is the possibility of raising possible explanations for such features.

3. 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 fractional 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 on the values of ϵ 0 , as we will show in the next subsections. The values of ρ , λ a , and λ b are assumed to remain constant during treatment. All rates (k, f, and h) are given in h 1 .
The values of the parameters were selected for simulating treatment conditions whose probability distributions governing the expression of pre-treated cells indicate qualitatively distinguishable steady state regimens as recently classified [23]. Figure 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 Figure 2, are computed by Equation (A33) and is denoted by ϕ ˜ n in the Appendix A.6.
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 quickly after 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 state duration is proportionally very long.
In next subsections, we present the results of simulations of the treatment response dynamics considering the five pre-treatment conditions shown in Figure 2. The trajectories of n ( t ) and σ ( t ) were obtained, respectively, using Equations (15) and (16) within the discrete intervals used to approximate the kinetic rates after treatment injection given by Equations (7) and (8).
Figure 3, Figure 4 and Figure 5 show the response after application of a single (or multiple) doses on graphs labeled as A–E (or F–J) followed by the number indicating the treatment scenario (1, 2, or 3). The pre-treatment conditions have ϵ 0 = ( 0.1, 1 , 2 , 10 , 10 ) . The black dashed lines indicate the aimed average value after treatment. The green lines indicate n ( t ) , and the red lines indicate the values of n ( t ) ± σ ( t ) . The parameters of the treatment are set to enable the average number of gene products to reach the post-treatment regime shown in Figure 2 for each set of parameter values.

3.1. Treatment Aiming at the OFF to ON Gene State Switching Rate

Figure 3 shows a simulation of the dynamics of the average number of mRNAs and its standard deviation resulting from a change on f after introduction of DETANONOate. The absolute error of each subinterval of the piece-wise approximation for f ( t ) is 1 × 10 4 . We assume the rates ( k , h , ρ ) are constant during treatment. For obtaining the Figure 3 A1–D1,F1–I1, we set ( k 0 , N 0 ) = ( k 1 , N 1 ) = ( 18.5, 110 ) and ( A 0 , A 1 ) = ( 0.09, 0.9) . Figure 3E1,J1 have ( k 0 , N 0 ) = ( k 1 , N 1 ) = ( 166.7, 1000 ) and ( A 0 , A 1 ) = ( 0.01, 0.1) to ensure the bursting gene expression regime during the pre-treatment stage.
On Figure 3A1–E1 (or Figure 3 F1–J1), we set f 0 = ( 0.0015, 0.015, 0.03, 0.15, 0.017) and f 1 = ( 0.14, 1.37, 2.73, 13.65, 0.18) , such that, h = ( 0.015, 0.15, 0.3, 1.52, 1.65) and ϵ 1 = ( 0.9, 9.1, 18.2, 91 , 11 ) with the given values of A 0 , A 1 and ρ . The maximum tolerated dose is considered to cause the steady state ON state probability to be multiplied by ten, as indicated by the values of A 0 and A 1 .
The first and second rows of graphs of Figure 3 show the simulation of the dynamics of expression after one and five drug doses, respectively, under five pre-treatment conditions.
Figure 3A1,B1 show the dynamics of response of the slow switching gene, a regime of synthesis of mRNAs whose numbers are governed by a bimodal distribution at the steady state. For ϵ 0 = 0.1 , we have the slowest treatment response and return to pre-treatment conditions.
The average number of mRNAs reaches a maximum of ∼80 after ∼30 h. The standard deviation is the second largest. Figure 3A1–D1 show the average mRNA number reaching or crossing the threshold in our simulations, the time for the average number to reach a maximum when ϵ 0 1 is ∼20 h, and the noise of the response decreases as we increase ϵ 0 . The exception is shown on Figure 3E1 where the noise is the largest and the time for the average number to reach a maximum is ∼10 h. In this case, the average number does not cross the threshold, and the lower line n ( t ) σ ( t ) does not appear.
The response to multi-dose treatment is shown in Figure 3F1–J1. The interval between doses was chosen to enable a sigmoidal-like response characterizing two distinct levels of expression with the average number increasing from 10 to at least ∼100. Figure 3F1 shows that the slow switching gene also causes the slowest response as n ( t ) crosses the threshold after ∼18 h. The curve for n ( t ) σ ( t ) also crosses the threshold after ∼35 h, which ensures a less-heterogeneous response to treatment. Figure 3G1–I1 show simulations for increasing the values of gene switching.
Those simulations show that n ( t ) crosses the threshold after ∼10 h, and the curves n ( t ) σ ( t ) reach higher values, which establishes the response to treatment with the minimal heterogeneity. Figure 3J1 shows the response of a burst gene, which n ( t ) crosses the threshold after ∼12 h and reaches a maximal value of ∼170. However, this regime causes the noisiest response as indicated by n ( t ) σ ( t ) not crossing the threshold. Note also that there are some “bumps” as the effect of the single dose is close to the maximum by the time of the next drug application. It is possible to prevent the bumps; however, we decided to show them to demonstrate that the time interval between drug dosages also requires one to consider the dynamics of the targeted gene.

3.2. Treatment Aiming at the RKIP mRNA Synthesis Rate

Figure 4 shows the dynamics of the average number of RKIP mRNAs and the standard deviations under treatment with 5-AzaC. The absolute error of each subinterval of the piece-wise approximation for k ( t ) is 1 × 10 4 . We assume that the rates ( f , h , ρ ) remain constant during treatment, which implies A 0 = A 1 and ϵ 0 = ϵ 1 also remain constant. For Figure 4A2–D2,F2–I2, we set ( k 0 , N 0 , k 1 , N 1 ) = ( 3.3, 20 , 33 , 200 ) , A 0 = 0.5 and f = h = ( 0.008, 0.08, 0.17, 0.83) . For Figure 4E2,J2, we set ( k 0 , N 0 , k 1 , N 1 ) = ( 166.7, 1000 , 1667 , 10 , 000 ) , A 0 = 0.01 and ( f , h ) = ( 0.017, 1.65) .
The first and second rows of graphs of Figure 4 show the simulation of the dynamics of expression after one and ten drug doses, respectively, under five pre-treatment conditions.
The treatments do not affect either ϵ or ρ , the decaying rates of the system. Then, Figure 4A2–E2 show that the average number of mRNAs reach the maximum (∼35 molecules) after similar intervals of ∼5 h and reach pre-treatment conditions after ∼40 h. Figure 4A2 shows the condition with the second largest standard deviation. Figure 4A2–D2 show that the noise of the response decreases with the increase of ϵ 0 . The burst regime shown in Figure 4E2 leads to the noisiest response. Indeed, n ( t ) σ ( t ) remains negative, while n ( t ) + σ ( t ) exceeds the threshold more than 2 × .
The response to multi-doses treatment is shown in Figure 4F2–J2. The interval between doses was chosen to enable a sigmoidal-like response with the average number reaching at least 100. All graphs show that n ( t ) crosses the threshold after ∼10 h . This is because the response to treatments and interval between drug application is the same in all scenarios. Since both ϵ 0 and ρ are not affected, the decaying rates to pre-treatment conditions are not affected.
In all scenarios, the curves for n ( t ) σ ( t ) do not cross the threshold because the treatment does not increase ϵ 0 that would cause a reduction in the noise in mRNA synthesis in a given scenario. However, Figure 4F2–I2 show that n ( t ) σ ( t ) reaches higher maximum values as ϵ 0 increases. Figure 4J2 shows the response in the case of transcriptional bursts, which leads to the most heterogeneous response to treatment as indicated by the curve n ( t ) + σ ( t ) reaching the highest maximum.

3.3. Treatment with the Two DRUGS Concomitantly

Figure 5 shows the dynamics of the average number of mRNAs and its standard deviation under treatment with both drugs: DETANONOate and 5-AzaC. The absolute errors of each subinterval of the piece-wise approximation for f ( t ) and k ( t ) are 1 × 10 4 for both. The rates ( h , ρ ) remain constant during treatment. For Figure 5A3–D3 (or Figure 5F3–I3), we set ( k 0 , N 0 , k 1 , N 1 ) = ( 18.5, 110 , 33 , 200 ) , ( A 0 , A 1 ) = ( 0.09, 0.5) , f 0 = ( 0.0015, 0.015, 0.03, 0.15) and f 1 = h = ( 0.015, 0.15, 0.3, 1.5) . For Figure 5E3,J3, we set ( k 0 , N 0 , k 1 , N 1 ) = ( 167 , 1000 , 333 , 2000 ) , ( A 0 , A 1 ) = ( 0.01, 0.05) , f 0 = 0.017 , f 1 = 0.087 and h = 1.65 . The values of A 0 , A 1 and ρ results in ϵ 1 = ( 0.18, 1.8, 3.6, 18 , 10.4) on Figure 5A3–E3 (or Figure 5F3–J3).
Figure 5 shows the simulations of the dynamics of response to treatment. Graphs of the first row indicate the results of simultaneous application of a single dose of both drugs. The second row shows the results of applications of multiple doses with drug targeting f (or k) being applied every 10 h (or 4 h) according to the agenda in Table 1.
Figure 5A3 shows the dynamics of the treatment response when the initial conditions are set for a slow switching gene. The average number reaches a maximum of ∼25 at t 35 h. The maximum of n ( t ) + σ ( t ) is ∼70, a highly heterogeneous response if we consider that the maximum standard deviation is equal to the maximum average—that is, a super-Poissonian regime of gene transcription. The return to pre-treatment is slow as indicated by the decay of the average after ∼200 h.
In Figure 5B3–D3, the maximum average numbers of mRNAs are ∼40, ∼45, and 50, respectively, reached after 15, 12, and 8 h. The curves n ( t ) σ ( t ) and n ( t ) + σ ( t ) become closer with the increase of ϵ 0 as shown from Figure 5A3–D3. The burst limit is different and, despite the high value of ϵ 0 , has larger noise as shown in Figure 5E3. At this limit, n ( t ) maximum is ∼45 reached after ∼8 h. n ( t ) + σ ( t ) reaches a maximum of ∼120, indicating the noisiest response. For all pre-treatment conditions, the average number does not cross the threshold.
The response to multi-doses treatment is shown in Figure 5F3–J3. Table 1 shows the sequences of applications of a given drug in fractions of maximal tolerated dose, ξ a and ξ b . The values were chosen aiming that the total amount of drugs is reduced in comparison with the single drug treatment. We also attempt to ensure a sigmoidal-like response such that the average number will reach at least ∼100. Figure 5F3 shows the dynamics of response for 8 (or 20) doses of drug targeting f (or k) with ξ a and ξ b ranging from 0.9 to 0.5 (details in Table 1).
The agenda enables a reduction of 20 % in comparison with application of full doses ( ξ a = ξ b = 1 ). As the pre-treatment gene is in a slow switching regime that leads to the slowest response as n ( t ) crosses the threshold after ∼40 h. It also has a high noise as can be noticed by the values of n ( t ) ± σ ( t ) . Figure 5G–J3 show the simulations of the response for 6 (or 15) doses of drugs a (or b) with ξ a and ξ b ranging from 0.8 to 0.5 (see the Table 1).
The cumulative reduction in full doses range from 29 % to 35 % (or from 23 % to 35 % ) in ξ a (or ξ b ). Figure 5G3–I3 show that n ( t ) crosses the threshold after ∼15 h. The curves n ( t ) ± σ ( t ) become closer as the value of ϵ increases, which indicates the direction to design treatment strategies with reduced response heterogeneity. Figure 5J3 shows the burst pre-treatment condition. Here, n ( t ) crosses the threshold after ∼22 h and n ( t ) + σ ( t ) reaches ∼200, which indicates the noisiest response.

3.4. Treatment Response Mapping ξ a and ξ b Fractional Effect of Drug Reduction

Figure 6 shows the treatment response measured in terms of the average (not weighted) values n σ ¯ (Figure 6A–E), and 2 σ ¯ (Figure 6F–J) from 60 to 80 h after the first dose. Five pre-treatment regimes labeled by ϵ 0 values are considered. The fraction of DETANONOate and 5-AzaC effect on f and k is indicated by ξ a and ξ b , respectively. We compute the average value of the function x ( t ) between instants t i and t f , denoted by x ( t ) ¯ , using the integral x ( t ) ¯ = 1 t f t i t i t f x ( t ) d t .
The absolute error of the piece-wise approximation for the exponential decay of the drug effect is set to 1 × 10 4 . To obtain each point of the heatmaps in Figure 6, we computed the dynamics of the treatment response to an agenda with doses with constant fractions ξ a and ξ b . DETANONOate was applied in eight doses with a time interval of 10 h, and 5-AzaC in 20 doses with an interval of 4 h. The grid is constructed by computing the trajectories for ξ a and ξ b varying from 0 to 1 by 0.05 increments.
Reduction of fractions of effects of dosage ξ a and ξ b for agendas shown in Table 1 and Figure 5 were based on heatmaps in Figure 6. The best scenario is located in the yellowest region of the heatmaps. The changes on the fractional effect of reduced doses over time are indicated by the sequence of red arrows within the graphs indicating each pre-treatment regime. Bimodal (Figure 6A,F) and burst (Figure 6E,J) regimes are the hardest pre-treatment conditions for dose reduction because there are almost no ( ξ a , ξ b ) ( 0 , 0 ) that ensure sufficiently high values of n σ ¯ (Figure 6A,E) and low 2 σ ¯ (Figure 6F,J).
The additional pre-treatment conditions, namely the remaining bimodal, the table-shaped, and the quasi-poisson (Figure 6B–D,G–I), enable better scenarios for dose reduction because of the larger values of ϵ 0 . Overall, all pre-treatment regimens enable a larger reduction of ξ a than of ξ b under the agenda proposed in Table 1 (see red arrows in Figure 6), because of the action of 5-AzaC on k (given in terms of ξ b ). Larger values of ξ b enable larger increases on average RKIP mRNA levels.
However, bimodal and burst regimens enable one to further reduce ξ b as the standard deviation remains constant or diminished, while the average RKIP mRNA numbers remain almost constant. However, though it is possible to reduce the variance, the average values do not reach sufficiently large values. Hence, we need to understand how treatment design may target the other kinetic rates, h and ρ , to enable an effective response independent of pre-treatment conditions.

3.5. Enhancing Ineffective Treatments Aiming at All Kinetic Rates

Response to enhanced treatment is shown in Figure 7 for five qualitatively different initial conditions aiming at optimal distribution. The hypothetical drug cocktail targets all rates of the model and we reach an optimal response time and heterogeneity reduction for all qualitatively distinguishable initial conditions, which, in previous subsections, led to unsatisfactory results. The pre-treatment distributions are shown by the red curves within the graphs in Figure 2: A2, B3, C2, D3, E3, which, respectively, describe a bimodal, bimodal limit, table shaped, quasi-Poissonian, and burst steady state regimes of gene transcription.
Our goal was to keep the average numbers of mRNAs around 100 and their fluctuations above the threshold. The decaying rates of the drugs affecting the drugs k, ρ , f, and h are denoted by λ i , where i indicates the rate targeted by the drug. The decaying rate of the drugs is fixed in h 1 : ( λ k , λ ρ , λ f , λ h ) = ( 0.25, 6 , 0.05, 0.053) , and the maximal tolerated dose of drug i is denoted by ξ i = 1 .
The values of λ f and λ k are the same used previously from DETANONOate and 5-AzaC, respectively, while λ ρ and λ h values were set based on the mRNA–microRNA binding duration time [88] and h in BACH1 half-life [89]. The time-dependence of h ( t ) and ρ ( t ) is described following the same framework used for k ( t ) and h ( t ) —see Equations (5)–(8). To obtain the dynamics of the average number of mRNAs and of the standard deviation, we extend the previous formulation of the master equation of Equations (10) and (11) to all kinetic rates.
The standard deviation of the post-treatment probability distribution is set to be
σ 1 = N 1 A 1 1 + N 1 ( 1 A 1 ) 1 + ϵ 1 ,
which has a local minimum at N 1 0 and a local maximum at
A 1 = 1 2 1 + 1 + ϵ 1 N 1 ,
where N is the maximal mRNA number, A is the steady state probability for the gene to be ON, ϵ is the gene switching speed (Equation (12)), and the subscript 1 indicates aimed post-treatment parameter values. A strategy to reduce σ 1 is to keep N 1 1 + ϵ 1 (or k 1 f 1 + h 1 + ρ 1 ), which ensures that σ 1 approaches the local maximum as A 1 1 . The reduction of σ 1 values depends on the mean number n 1 desired, so that n 1 A 1 ( 1 + ϵ 1 ) .
To enhance the previously ineffective treatments designs, we used the aforementioned optimization approach. The drug dose and interval of application are set to enable ϵ 1 = 91 , which was found to be the value leading to the successful treatment design shown in D1 of Figure 2. Indeed, this resulted in the best response dynamics of the average number of RKIP mRNA’s as shown in I1 of Figure 3.
Figure 7A4 and B4 show the dynamics of response to enhanced treatment when pre-treatment RKIP gene transcription is governed by a slow switching promoter giving rise to bimodal probability distributions ( ϵ 0 1 ). The previously ineffective treatment design aimed at bimodal and table-shaped regimens are, respectively, shown by the green lines of A2 and B3 in Figure 2. The dynamics of response to the enhanced treatment designs are shown in Figure 7A2*,B3*. Both optimized responses resulted from similar changes in kinetic parameters that caused an increase in A 1 and a reduction in N 1 and have similar agendas (detailed in Appendix B).
Enhanced treatment design aiming for a pre-treatment fast switching gene shown in C4 and D4 of Figure 7. The initial probability distributions are table-shaped and quasi-Poissonian and have ϵ 0 > 1 . The response dynamics to ineffective treatment designs are shown in C2 and D3 of Figure 2. The response to enhanced treatment with reduced heterogeneity and increased speed is shown in C2* and D3* of Figure 7. The enhanced treatment agendas were the same in both cases (see parameters and agendas in Appendix B). The value of n 1 was reduced in comparison to the ineffective treatment design because of the decrease of N 1 not being compensated by the increase of A 1 .
The ineffective treatment design for an initially burst regimen of RKIP transcription shown in Figure 2E3 is enhanced. The resulting post-treatment distribution is a quasi-Poissonian. The dynamics of response to enhanced treatment design is shown, respectively, in E3* and E4 of Figure 7. This enhancement required the largest reductions in n 1 in comparison to the ineffective design by increasing A 1 and decreasing N 1 . Note, however, the clear reduction in the standard deviations (heterogeneity) of the response. The parameters values and agenda applied are detailed in Appendix B.

4. Discussion

When the binary model for regulation of gene expression is considered, it is useful to characterize the pre-treatment gene expression regime in terms of A 0 , ϵ 0 , and N 0 (see Equation (12)). Those rates will affect the pre-existing (and post-treatment) noise on the number of mRNAs, see Equations (13) and (14) (and Figure 3, Figure 4 and Figure 5 and Figure 7). Furthermore, the time of response to treatment also depends on ϵ 0 , because the decaying to steady state regime of the average number of mRNAs and its variance have both ρ and ϵ ρ as their smallest rates (see Equation (15)). For ϵ > 1 and < 1 the time for the system to reach the steady state is, respectively, ρ 1 and ϵ ρ 1 .
All trajectories of the average values shown in graphs of Figure 4 reach a maximal (Figure 4A2–E2) or cross the threshold (Figure 4F2–J2) after an approximately fixed interval. This is because the component e ϵ ρ t of n ( t ) is null as the treatment does not change the steady state probability for the gene to be ON (see Equation (A20)). Additionally, the noise on treatment response, measured in terms of the trajectories of n ( t ) ± σ ( t ) , is larger when we have the smaller values of the relative switching speed. Inspection of trajectories shown in A to D and F to I of Figure 3, Figure 4 and Figure 5 helps to identify those features.
The pre-treatment condition, which enables the fastest and least heterogeneous response takes place for higher values of the relative switching speed and for the probability for the gene to be ON being >0.1. D and I of Figure 3, Figure 4 and Figure 5 show that those are the conditions enabling the smallest differences between trajectories n ( t ) ± σ ( t ) . At this limit, the gene is behaving as a quasi-Poissonian source of transcripts [23]. As the value of the pre-treatment switching speed is reduced from toward two, we have an increase in the noise of the treatment response, though the speed of response does not increase significantly, as shown in C and H of Figure 3, Figure 4 and Figure 5.
The burst pre-treatment regime ( ϵ 0 = 10 ) leads to the noisiest responses to treatment, although being among the fastest ones. The noisy response to treatment occurs partially because when we simulated treatment scenarios shown in E and J of Figure 3, Figure 4 and Figure 5, the value of the A ( t ) remained small. In that case, the dynamics would still remain in a transcriptional burst regime. The small increase in A ( t ) is because we assumed that the maximal tolerated dose of the drug responsible for increasing f would cause the probability of the ON state to be 10× larger, independent of its pre-treatment value. Such an assumption needs to be confirmed by experimental studies, and an alternative formulation might be proposed because of a lack of confirmation. B to E and G to J of Figure 3, Figure 4 and Figure 5 show that the time of response to treatment are all similar since the rates of decay to steady state are ≤ ρ for ϵ 1 .
The slow relative switching speed pre-treatment condition (here, we set ϵ 0 = 0.1, 1 ) establishes a challenging regime to be approached. The response is slower and noisier if one takes the strategies considered here. A and F of Figure 3 and Figure 5 indicate that the increase in the average number of mRNA is the slowest when we have the minimal ϵ despite treatment increasing f will cause an increase of the relative switching speed.
In the high noise responses, the trajectories n ( t ) σ ( t ) do not cross the threshold. This indicates that, in a population of cells with the same regulatory conditions of expression of RKIP gene, the response to treatment is heterogeneous and may be insufficient. Therefore, the dosage and target need to be constructed to ensure a reduction in the response heterogeneity. This can be carried out by increasing the switching speed such that the gene expression will approach a quasi-Poissonian regime. Increasing the switching speed also has the benefit of reducing the time of the response to therapy.
The fractional effect of dose reduction of DETATANOate and 5-AzaC shown in Figure 6 indicates the difficulties of optimizing the agenda for pre-treatment regimes far from quasi-Poissonian. We highlight bimodal and burst regimes, where the responses need to have larger increases regarding the average mRNA levels. For that, one needs to target all effective kinetic rates of the stochastic binary model for gene expression. For the case of RKIP, one might also target Snail, BACH1, Sp1, CREB, and p300 repressors to change h, while a variety of miRNAs and their possible inactivators would modulate ρ [32].
The responses to enhanced treatments shown in Figure 7 (Figure 7A4–E4) indicate the potential of our approach to help on the design of low dose multi-drugs administered in pulses. The three qualitatively different pre-treatment regimens, namely slow bimodal, fast quasi-Poissonian, and burst, have distinguishable response curves. Those are resulting from the specificities of the enhanced treatment design aiming at keeping the average mRNA numbers and deviations (heterogeneity) above the threshold.
The synthesis of RKIP protein may have a non-linear dependence on the mRNA numbers because translation can be regulated by mechanisms, such as mRNA stability, translational control, and proteosomal degradation [78]. Hence, a theoretical approach considering translation might require the addition of two or more characteristic timescales depending on protein synthesis and degradation [90]. We emphasize that the amounts of RKIP mRNAs in our model are the net effect of the enhanced treatment, since its expression results from the interaction with a complex gene network. Our approach, however, provides the building blocks of gene networks [19] that can be used to understand the functioning of larger modules [14,15,16] and, hence, to further enhance treatment designs of cancer at metastatic stage.
5-AzaC and DETANONOate are non-specific drugs to increase RKIP, which may also interfere with additional biological processes taking place in tumoral or normal cells. The current study does not provide an insight on those problems and specific experimental designs are needed to approach these. Initially, one may consider simulations using pre-treatment conditions corresponding to the normal condition. In that case, it is fair to consider that the fraction of effect of reduced dose is smaller, if we assume that the drugs have a higher absorption in cancer cells. A more complex model for the dynamics of multiple genes is needed to investigate the overall effect of non-specific treatment designs. Alternatively, one might design treatments with higher specificity to further investigate the validity of our approach for simulating gene therapeutic designs.
Another possible application of our approach is on the investigation of therapeutic designs based on 5-AzaC and DETANONOate. For example, 5-AzaC promotion of DNA demethylation could be used to regulate the expression of genes involved in reverting adult stem cell ageing [91], while its inhibitory properties of DNA methyltransferase has the potential for regenerating mature mammalian inner ear hair cells [92]. NO donors as DETANONOate are involved in many cell processes, such as vasodilation, neurotransmission, macrophage-mediated immunity, and anti-inflammatory responses [93]. One application of our approach is to help in elucidating the mechanisms of those therapeutic designs.
To further improve our approach, one needs to perform a data-based validation of our model. One possible experiment would involve the decrease of RKIP levels using a siRNA approach to simulate different levels of activation of the RKIP gene. This would be a more specific and precise way to interfere in the RKIP pathway compared with using 5-AzaC. We could make different levels of RKIP inhibition by siRNA and see if the modeling can predict MEK and ERK phosphorylation levels, the expression of MAPK-induced genes, and the consequence of these activations in cell migration, proliferation, and invasion. The siRNA approach is not able to check the OFF state of the RKIP gene; however, it is possible to use a CRISPR/CAS9 system to eliminate the RKIP gene, and then have the OFF state information. Inference methods would be employed on the estimation of parameter values and RKIP mRNA numbers.
Those experiments would help to understand how environmental information is processed by cells by investigating how regulation of expression of RKIP gene affects its relative downstream kinases. In recent work, we demonstrated that the slow switching speed regime maximizes the mutual information between the number of transcripts and promoter activity [34]. This regime also coincides with conditions of highly heterogeneous and slow response to treatment. Hence, we may ask whether it is relevant for a master regulatory gene to transmit information about its own promoter state and in which cellular context, namely, in a cancerous or healthy one. The answer to this question is important because it helps us to understand the role of noise and the conditions under which its reduction is desirable [33] as it happens when we have a negatively self-regulating gene.

5. Conclusions

In this manuscript, we present a stochastic binary model for transcription of the RKIP gene with treatment-induced time-dependent kinetic rates. The exact solutions of this model are approximated using the exact solutions of the equivalent model with stationary kinetic rates. This is the simplest exactly solvable model to describe the regulated transcription of the RKIP gene. This enables us to simulate the effects of the application of a drug that changes one (or multiple) kinetic parameters participating in the regulation of RKIP gene.
To demonstrate the usefulness of our approach, we simulated three scenarios in which we aimed to increase the number of RKIP mRNAs by increasing the: i. OFF to ON switching rate using DETANONOate; ii. synthesis rate using 5-AzaC; and iii. both rates of a gene using both drugs together. We showed that treatment response speed and heterogeneity depended on the pre-treatment state of the gene. Then, we presented an enhanced treatment design that ensured reduced heterogeneity and time of response. In addition to being useful to inspect treatment designs, the response to treatment may be used for inference of the kinetic constants of a given gene in a synthetic system.

Author Contributions

Conceptualization, A.F.R.; methodology, A.F.R.; software, G.G., L.R.G. and A.F.R.; validation, G.G., L.R.C.B., L.R.G., T.C.T.J. and A.F.R.; formal analysis, G.G., L.R.C.B., L.R.G., T.C.T.J. and A.F.R.; investigation, G.G. and A.F.R.; resources, A.F.R.; data curation, G.G., L.R.G. and A.F.R.; writing—original draft preparation, A.F.R.; writing—review and editing, G.G., L.R.C.B., L.R.G., T.C.T.J. and A.F.R.; visualization, G.G. and A.F.R.; supervision, A.F.R.; project administration, A.F.R.; funding acquisition, A.F.R. All authors have read and agreed to the published version of the manuscript.

Funding

The authors thank funders for support: G.G. thanks CAPES and PPG-Modelagem de Sistemas Complexos (88887.506433/2020-00); L.R.C.B. thanks PRONON Retratos da Mama (25000.069252/2015-79); T.C.T.J. thanks FAPESP (2015/22814-5); L.R.G. thanks Programa Unificado de Bolsas (PUB-USP) and Pró-Reitoria de Graduação da USP; and A.F.R. thanks CAPES (88881.062174/2014-01), FAPESP (2012/24962-3), and Magalu Initiative for Global Health.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The codes can be provided under request to the corresponding author.

Acknowledgments

The autors would like to thank Marsha R. Rosner and Roger Chammas for stimulating and helpful discussions that contributed to the concept behind this story. The autors thank anonymous reviewers for comments which enhanced the quality of our manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Formulae

Appendix A.1. Piecewise Approximation for Time-Dependent Kinetic Parameters

We denote by E i the error of the piece-wise approximation for the decay of the drug effect on a given rate within an interval [ t i , t i + δ i ) . We compute the error as the absolute value of the difference between area of the rectangle with base δ i > 0 and height e λ t i , and the area under the exponential function e λ t within the interval [ t i , t i + δ i ] , namely:
E i = δ i e λ t i t i t i + δ i e λ t d t = e λ t i λ λ δ i 1 + e λ δ i .
Note that E i is positive definite.
Once the error E i on an interval is set, we may solve a transcendental equation λ δ i + 1 + λ e λ t i E i = e λ δ i for δ i using some numerical procedure. In this manuscript, we set E i = 10 4 for all subintervals along [ 0 , t max ] .

Appendix A.2. System of Coupled ODEs Governing the Moments of the Distribution

Here, we describe the process for obtaining the time-dependent solutions for n ( t ) and σ 2 ( t ) given at Equations (15) and (16). This is carried out considering the kinetic rates to be constant in the master Equations (10) and (11), which become:
d α n d t = k ( α n 1 α n ) + ρ [ ( n + 1 ) α n + 1 n α n ] h α n + f β n ,
d β n d t = ρ [ ( n + 1 ) β n + 1 n β n ] + h α n f β n .
The master equations above govern the dynamics of the joint probability distribution Π ( α n , β n ) . This distribution has two marginal probability distributions: i. the probability of finding n mRNAs, independently of the promoter state, is denoted by ϕ n and given by:
ϕ n ( t ) = α n ( t ) + β n ( t ) ;
ii. The probability of finding the promoter state as ON or OFF is, respectively, denoted by A ( t ) and B ( t ) and given by:
A ( t ) = n = 0 α n , and B ( t ) = n = 0 β n .
In what follows, we describe how the generating functions technique can be used to obtain explicit expressions for these probabilities and associated moments.
The master Equations (A2) and (A3) can be viewed as an infinite system of coupled ordinary differential equations whose exact solution is usually difficult to obtain. In solving it, we may employ the generating functions technique, which considers the probabilities α n ( t ) and β n ( t ) as coefficients of a Taylor expansion of an analytic function, namely:
α ( z , t ) = n = 0 α n ( t ) z n β ( z , t ) = n = 0 β n ( t ) z n .
The application of the definitions above enables one to transform Equations (A2) and (A3) into a pair of partial differential equations (PDEs) on variables ( z , t ) , namely:
α t = ( z 1 ) k α ρ α z h α + f β ,
β t = ρ ( z 1 ) β z + h α f β .
This pair of coupled PDEs is obtained by first multiplying both sides of Equations (A2) and (A3) by z n , and then, by application of the operator n = 0 .
Here, we are interested in understanding the dynamics of both n ( t ) and σ 2 ( t ) . For that, we need to recover the system of ODEs governing the p-th order moments of the distribution Π ( α n , β n ) , where p = 1 , 2 , . This is done by applying the z z p on Equations (A7) and (A8) and, on the resulting expression, setting z = 1 . This procedure is proposed because the p-th order moment may be recovered from the generating functions as follows:
n α p = n = 0 α n n p = z z p α ( z , t ) z = 1 ,
n β p = n = 0 β n n p = z z p β ( z , t ) z = 1 .
n p = n = 0 ϕ n n p = z z p ϕ ( z , t ) z = 1 .
Here, ϕ ( z , t ) = α ( z , t ) + β ( z , t ) is the generating function of the marginal probability of finding n mRNAs, ϕ n . Hence, n p is the p-th order moment on the number of mRNAs.
The expressions of n ( t ) and σ 2 ( t ) are obtained after solving the following system composed of linear non-homogeneous coupled ODEs:
d A d t + f + h A = f ,
d n d t + ρ n = k A ,
d n α d t + f + h + ρ n α = k A + f n ,
d n 2 d t + 2 ρ n 2 = k A + ρ n + 2 k n α ,
where we omitted the time-dependence of A, n , n α , and n 2 .
The above system of ODEs was obtained as follows. 1. Equation (A12): at Equation (A7) we set, in the following order, β = ϕ α , z = 1 , α 1 , t = A t , and ϕ 1 , t = 1 . 2. Equation (A13): in following order, we sum Equations (A7) and (A8), set β = ϕ α , apply the operator z , evaluate the resulting equation at z = 1 , and use the identities α 1 , t = A t , and ϕ z z = 1 = n . 3. Equation (A14): at Equation (A7) we set, in the following order, β = ϕ α , apply the operator z , evaluate the resulting equation at z = 1 , and use the identities α 1 , t = A t , ϕ z z = 1 = n , and α z z = 1 = n α . 4. Equation (A15): in following order, we sum Equations (A7) and (A8), set β = ϕ α , apply the operator 2 z 2 , evaluate the resulting equation at z = 1 , and use the identities α 1 , t = A t , ϕ z z = 1 = n , α z z = 1 = n α , and 2 ϕ z 2 z = 1 = n 2 n .

Appendix A.3. Steady State Values of the Moments of the Distribution

To obtain the steady state solutions of the system of EDOs of Equations (A12)–(A15) one may set A A s , n n s , n α n α s , and n 2 n 2 s , where the index s indicates these quantities at the steady state limit, when they become constant. Hence, the derivatives are null and one may solve the resulting system of equations. Equation (A12) is solved for A s , which is replacedinthe next, and so on. This procedure will result in
A s = f f + h ,
n s = A s N
n α s = A s N 1 + A s ϵ 1 + ϵ ,
n 2 s = A s N 1 + N 1 + A s ϵ 1 + ϵ .

Appendix A.4. Coefficients of Equations (15) and (16)

Y = N A 0 A s 1 ϵ ,
V = n 0 Y n s ,
U 1 = U 2 Y n s , W 1 = W 2 Y V , X 1 = X V 2 ,
U = Y 1 + 2 N 1 ϵ 1 A s 2 ϵ ,
W = 2 N 1 ϵ N 1 A s ϵ A s A 0 A 0 ϵ + 1 A s n 0 + n α 0 ,
X = n 2 0 n 2 s U 1 + 2 n s V W .

Appendix A.5. Exact Formulas for A(t), 〈nα〉(t), and 〈n2〉(t)

The exact time-dependent solutions of the system of EDOs of Equations (A12)–(A15) can be obtained using the method of variation of the constant (see Section 9.2 of [94]). For that, one may solve Equation (A12) for A ( t ) , insert it in the next equation, and so on. The resulting solutions are:
A ( t ) = A s + A 0 A s e h + f t ,
n ( t ) = n s + Y e ϵ ρ t + V e ρ t ,
n α ( t ) = n α s + R e h + f t + S e ρ t + T e f + h + ρ t ,
n 2 ( t ) = n 2 s + U e h + f t + 1 + 2 n s V e ρ t + W e f + h + ρ t + X e 2 ρ t ,
where
R = Y 1 ϵ 1 A s ,
S = A s n 0 + n s A s ϵ A 0 1 ϵ ,
T = n α 0 n α s R S .

Appendix A.6. Steady State Probabilities of Finding n Gene Products

The steady state probability distributions presented in this manuscript were obtained from the exact solutions of the stochastic binary model for regulation of gene expression presented in [19,95]. Here, we denote the steady state marginal probabilities of finding n gene products by ϕ ˜ n . It is defined in terms of the probabilities of finding the gene ON or OFF when n gene products are found inside the cell, which are denoted by α ˜ n or β ˜ n . Hence we write ϕ ˜ n = α ˜ n + β ˜ n , which is given by:
ϕ ˜ n ( A s , ϵ , N ) = N n n ! ( A s ϵ ) n ( ϵ ) n M ( A s ϵ + n , ϵ + n , N ) ,
where ( a ) n a ( a + 1 ) ( a + n 1 ) is the Pochhammer symbol and M ( a , b , z ) = m = 0 ( a ) m ( b ) m z m m ! is the KummerM function ([96], p. 503).

Appendix B. Treatment Parameters and the Agenda of Doses for Enhanced Treatment Designs

The following table shows the treatment parameters used to obtain the graphs of Figure 7. The post-treatment kinetic rates, ( k 1 , ρ 1 , f 1 , h 1 ) , and distributions parameters, ( N 1 , A 1 , n 1 , σ 1 ) , where set to enable the enhanced treatment designs. The number of doses is represented by n i and the time interval between them by τ i , where i indicates the aimed kinetic rate. For all treatments, ϵ 1 = 91 , n f = 3 and τ f = 10 .
Graph ( k 1 , ρ 1 , f 1 , h 1 ) ( N 1 , A 1 , n 1 , σ 1 ) ( n k , n ρ , n h ) ( τ k , τ ρ , τ h )
A4 ( 112 , 0.83, 58 , 18 ) ( 135 , 0.76, 103 , 11.8) ( 8 , 30 , 30 ) ( 4 , 1 , 1 )
B4 ( 91 , 0.67, 40 , 21 ) ( 136 , 0.65, 88 , 11.6) ( 8 , 30 , 15 ) ( 4 , 1 , 2 )
C4/D4 ( 44 , 0.42, 23 , 15 ) ( 106 , 0.6, 63.6, 9.6) ( 30 , 30 , 30 ) ( 1 , 1 , 1 )
E4 ( 235 , 1.17, 23 , 83 ) ( 201 , 0.22, 44 , 10.9) ( 8 , 10 , 10 ) ( 4 , 3 , 3 )

References

  1. Bulaklak, K.; Gersbach, C.A. The once and future gene therapy. Nat. Commun. 2020, 11, 5820. [Google Scholar] [CrossRef] [PubMed]
  2. Thakore, P.I.; Black, J.B.; Hilton, I.B.; Gersbach, C.A. Editing the epigenome: Technologies for programmable transcription and epigenetic modulation. Nat. Methods 2016, 13, 127–137. [Google Scholar] [CrossRef]
  3. Cai, W.; Zhou, W.; Han, Z.; Lei, J.; Zhuang, J.; Zhu, P.; Wu, X.; Yuan, W. Master regulator genes and their impact on major diseases. PeerJ 2020, 8, e9952. [Google Scholar] [CrossRef] [PubMed]
  4. Yesilkanal, A.E.; Yang, D.; Valdespino, A.; Tiwari, P.; Sabino, A.U.; Nguyen, L.C.; Lee, J.; Xie, X.H.; Sun, S.; Dann, C.; et al. Limited inhibition of multiple nodes in a driver network blocks metastasis. eLife 2021, 10, e59696. [Google Scholar] [CrossRef]
  5. Yesilkanal, A.E.; Johnson, G.L.; Ramos, A.F.; Rosner, M.R. New strategies for targeting kinase networks in cancer. J. Biol. Chem. 2021, 297, 101128. [Google Scholar] [CrossRef] [PubMed]
  6. Delbrück, M. Statistical Fluctuations in Autocatalytic Reactions. J. Chem. Phys. 1940, 8, 120–124. [Google Scholar] [CrossRef]
  7. Marusyk, A.; Almendro, V.; Polyak, K. Intra-tumour heterogeneity: A looking glass for cancer? Nat. Rev. Cancer 2012, 12, 323–334. [Google Scholar] [CrossRef]
  8. Michor, F.; Beal, K. Improving Cancer Treatment via Mathematical Modeling: Surmounting the Challenges Is Worth the Effort. Cell 2015, 163, 1059–1063. [Google Scholar] [CrossRef] [Green Version]
  9. Alizadeh, A.A.; Aranda, V.; Bardelli, A.; Blanpain, C.; Bock, C.; Borowski, C.; Caldas, C.; Califano, A.; Doherty, M.; Elsner, M.; et al. Toward understanding and exploiting tumor heterogeneity. Nat. Med. 2015, 21, 846–853. [Google Scholar] [CrossRef]
  10. Brock, A.; Krause, S.; Ingber, D.E. Control of cancer formation by intrinsic genetic noise and microenvironmental cues. Nat. Rev. Cancer 2015, 15, 499–509. [Google Scholar] [CrossRef]
  11. Welch, D.R. Tumor Heterogeneity—A ‘Contemporary Concept’ Founded on Historical Insights and Predictions. Cancer Res. 2016, 76, 4–6. [Google Scholar] [CrossRef] [Green Version]
  12. Guinn, M.T.; Wan, Y.; Levovitz, S.; Yang, D.; Rosner, M.R.; Balázsi, G. Observation and Control of Gene Expression Noise: Barrier Crossing Analogies Between Drug Resistance and Metastasis. Front. Genet. 2020, 11, 586726. [Google Scholar] [CrossRef] [PubMed]
  13. Bhalla, U.S.; Iyengar, R. Emergent Properties of Networks of Biological Signaling Pathways. Science 1999, 283, 381–387. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  14. Balázsi, G.; van Oudenaarden, A.; Collins, J.J. Cellular decision making and biological noise: From microbes to mammals. Cell 2011, 144, 910–925. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  15. Ferrell, J.E. Perfect and Near-Perfect Adaptation in Cell Signaling. Cell Syst. 2016, 2, 62–67. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  16. Gérard, C.; Gonze, D.; Goldbeter, A. Revisiting a skeleton model for the mammalian cell cycle: From bistability to Cdk oscillations and cellular heterogeneity. J. Theor. Biol. 2019, 461, 276–290. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  17. Purvis, J.E.; Lahav, G. Encoding and Decoding Cellular Information through Signaling Dynamics. Cell 2013, 152, 945–956. [Google Scholar] [CrossRef] [Green Version]
  18. Levine, J.H.; Lin, Y.; Elowitz, M.B. Functional Roles of Pulsing in Genetic Circuits. Science 2013, 342, 1193–1200. [Google Scholar] [CrossRef] [Green Version]
  19. Ramos, A.F.; Innocentini, G.C.P.; Forger, F.M.; Hornos, J.E.M. Symmetry in biology: From genetic code to stochastic gene regulation. IET Syst. Biol. 2010, 4, 311–329. [Google Scholar] [CrossRef] [Green Version]
  20. Peccoud, J.; Ycart, B. Markovian modelling of gene product synthesis. Theor. Popul. Biol. 1995, 48, 222–234. [Google Scholar] [CrossRef]
  21. Iyer-Biswas, S.; Hayot, F.; Jayaprakash, C. Stochasticity of gene products from transcriptional pulsing. Phys. Rev. E 2009, 79, 031911. [Google Scholar] [CrossRef] [Green Version]
  22. Ramos, A.F.; Gama, L.R.; Morais, M.C.C.; Martins, P.C.M. Stochastic modeling for investigation of the regulation of transcription of the RKIP gene. In Prognostic and Therapeutic Applications of RKIP in Cancer; Bonavida, B., Baritaki, S., Eds.; Elsevier: Amsterdam, The Netherlands, 2020; pp. 257–276. [Google Scholar] [CrossRef]
  23. Giovanini, G.; Sabino, A.U.; Barros, L.R.C.; Ramos, A.F. A comparative analysis of noise properties of stochastic binary models for a self-repressing and for an externally regulating gene. Math. Biosci. Eng. 2020, 17, 5477–5503. [Google Scholar] [CrossRef]
  24. Jarrett, A.M.; Faghihi, D.; Hormuth, D.A.; Lima, E.A.B.F.; Virostko, J.; Biros, G.; Patt, D.; Yankeelov, T.E. Optimal Control Theory for Personalized Therapeutic Regimens in Oncology: Background, History, Challenges, and Opportunities. J. Clin. Med. 2020, 9, 1314. [Google Scholar] [CrossRef]
  25. Bonavida, B.; Baritaki, S. (Eds.) Prognostic and Therapeutic Applications of RKIP in Cancer; Elsevier: Amsterdam, The Netherlands, 2020. [Google Scholar] [CrossRef]
  26. Martinho, O.; Granja, S.; Jaraquemada, T.; Caeiro, C.; Miranda-Gonçalves, V.; Honavar, M.; Costa, P.; Damasceno, M.; Rosner, M.R.; Lopes, J.M.; et al. Downregulation of RKIP Is Associated with Poor Outcome and Malignant Progression in Gliomas. PLoS ONE 2012, 7, e30769. [Google Scholar] [CrossRef] [Green Version]
  27. Martinho, O.; Pinto, F.; Granja, S.; Miranda-Gonçalves, V.; Moreira, M.A.R.; Ribeiro, L.F.J.; di Loreto, C.; Rosner, M.R.; Longatto-Filho, A.; Reis, R.M. RKIP Inhibition in Cervical Cancer Is Associated with Higher Tumor Aggressive Behavior and Resistance to Cisplatin Therapy. PLoS ONE 2013, 8, e59104. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  28. Lamiman, K.; Keller, J.M.; Mizokami, A.; Zhang, J.; Keller, E.T. Survey of Raf Kinase Inhibitor Protein (RKIP) in Multiple Cancer Types. Crit. Rev. Oncog. 2014, 19, 455–468. [Google Scholar] [CrossRef] [PubMed]
  29. Yesilkanal, A.; Rosner, M. Targeting Raf Kinase Inhibitory Protein Regulation and Function. Cancers 2018, 10, 306. [Google Scholar] [CrossRef] [Green Version]
  30. Dangi-Garimella, S.; Yun, J.; Eves, E.M.; Newman, M.; Erkeland, S.J.; Hammond, S.M.; Minn, A.J.; Rosner, M.R. Raf kinase inhibitory protein suppresses a metastasis signalling cascade involving LIN28 and let-7. EMBO J. 2009, 28, 347–358. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  31. Fu, Z.; Smith, P.C.; Zhang, L.; Rubin, M.A.; Dunn, R.L.; Yao, Z.; Keller, E.T. Effects of Raf Kinase Inhibitor Protein Expression on Suppression of Prostate Cancer Metastasis. J. Natl. Cancer Inst. 2003, 95, 878–889. [Google Scholar] [CrossRef]
  32. Zaravinos, A.; Bonavida, B.; Chatzaki, E.; Baritaki, S. RKIP: A Key Regulator in Tumor Metastasis Initiation and Resistance to Apoptosis: Therapeutic Targeting and Impact. Cancers 2018, 10, 287. [Google Scholar] [CrossRef] [Green Version]
  33. Ramos, A.F.; Reinitz, J. Physical implications of so(2, 1) symmetry in exact solutions for a self-repressing gene. J. Chem. Phys. 2019, 151, 041101. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Gama, L.R.; Giovanini, G.; Balázsi, G.; Ramos, A.F. Binary Expression Enhances Reliability of Messaging in Gene Networks. Entropy 2020, 22, 479. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  35. Yeung, K.; Seitz, T.; Li, S.; Janosch, P.; McFerran, B.; Kaiser, C.; Fee, F.; Katsanakis, K.D.; Rose, D.W.; Mischak, H.; et al. Suppression of Raf-1 kinase activity and MAP kinase signalling by RKIP. Nature 1999, 401, 173–177. [Google Scholar] [CrossRef] [Green Version]
  36. Trakul, N.; Menard, R.E.; Schade, G.R.; Qian, Z.; Rosner, M.R. Raf Kinase Inhibitory Protein Regulates Raf-1 but Not B-Raf Kinase Activation. J. Biol. Chem. 2005, 280, 24931–24940. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  37. Yesilkanal, A.E.; Rosner, M.R. Raf Kinase Inhibitory Protein (RKIP) as a Metastasis Suppressor: Regulation of Signaling Networks in Cancer. Crit. Rev. Oncog. 2014, 19, 447–454. [Google Scholar] [CrossRef] [Green Version]
  38. Yeung, K.C.; Rose, D.W.; Dhillon, A.S.; Yaros, D.; Gustafsson, M.; Chatterjee, D.; McFerran, B.; Wyche, J.; Kolch, W.; Sedivy, J.M. Raf Kinase Inhibitor Protein Interacts with NF-κB-Inducing Kinase and TAK1 and Inhibits NF-κB Activation. Mol. Cell. Biol. 2001, 21, 7207–7217. [Google Scholar] [CrossRef] [Green Version]
  39. Zhao, J.; Wenzel, S. Interactions of RKIP with Inflammatory Signaling Pathways. Crit. Rev. Oncog. 2014, 19, 497–504. [Google Scholar] [CrossRef] [Green Version]
  40. Lorenz, K.; Lohse, M.J.; Quitterer, U. Protein kinase C switches the Raf kinase inhibitor from Raf-1 to GRK-2. Nature 2003, 426, 574–579. [Google Scholar] [CrossRef]
  41. Corbit, K.C.; Trakul, N.; Eves, E.M.; Diaz, B.; Marshall, M.; Rosner, M.R. Activation of Raf-1 Signaling by Protein Kinase C through a Mechanism Involving Raf Kinase Inhibitory Protein. J. Biol. Chem. 2003, 278, 13061–13068. [Google Scholar] [CrossRef] [Green Version]
  42. Shvartsur, A.; Givechian, K.B.; Garban, H.; Bonavida, B. Overexpression of RKIP and its cross-talk with several regulatory gene products in multiple myeloma. J. Exp. Clin. Cancer Res. 2017, 36, 62. [Google Scholar] [CrossRef] [Green Version]
  43. Datar, I.; Tegegne, H.; Qin, K.; Al-Mulla, F.; Bitar, M.S.; Trumbly, R.J.; Yeung, K.C. Genetic and Epigenetic Control of RKIP Transcription. Crit. Rev. Oncog. 2014, 19, 417–430. [Google Scholar] [CrossRef] [PubMed]
  44. Galal, Y.; Zaravinos, A.; Bonavida, B. Regulation of NKG2D by RKIP: Implications on NK-mediated cytotoxicity and cytokine production. In Successes and Challenges of NK Immunotherapy; Elsevier: Amsterdam, The Netherlands, 2021; pp. 233–265. [Google Scholar] [CrossRef]
  45. Li, D.X.; Cai, H.Y.; Wang, X.; Feng, Y.L.; Cai, S.W. Promoter methylation of Raf kinase inhibitory protein: A significant prognostic indicator for patients with gastric adenocarcinoma. Exp. Ther. Med. 2014, 8, 844–850. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  46. Wei, H.; Liu, Z.; She, H.; Liu, B.; Gu, J.; Wei, D.; Zhang, X.; Wang, J.; Qi, S.; Ping, F. Promoter methylation and expression of Raf kinase inhibitory protein in esophageal squamous cell carcinoma. Oncol. Lett. 2017, 13, 1866–1872. [Google Scholar] [CrossRef] [Green Version]
  47. Beach, S.; Tang, H.; Park, S.; Dhillon, A.S.; Keller, E.T.; Kolch, W.; Yeung, K.C. Snail is a repressor of RKIP transcription in metastatic prostate cancer cells. Oncogene 2007, 27, 2243–2248. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  48. Labbozzetta, M.; Poma, P.; Vivona, N.; Gulino, A.; D’Alessandro, N.; Notarbartolo, M. Epigenetic changes and nuclear factor-κB activation, but not microRNA-224, downregulate Raf-1 kinase inhibitor protein in triple-negative breast cancer SUM 159 cells. Oncol. Lett. 2015, 10, 3807–3815. [Google Scholar] [CrossRef] [Green Version]
  49. Ren, G.; Baritaki, S.; Marathe, H.; Feng, J.; Park, S.; Beach, S.; Bazeley, P.S.; Beshir, A.B.; Fenteany, G.; Mehra, R.; et al. Polycomb Protein EZH2 Regulates Tumor Invasion via the Transcriptional Repression of the Metastasis Suppressor RKIP in Breast and Prostate Cancer. Cancer Res. 2012, 72, 3091–3104. [Google Scholar] [CrossRef] [Green Version]
  50. Lee, J.; Lee, J.; Farquhar, K.S.; Yun, J.; Frankenberger, C.A.; Bevilacqua, E.; Yeung, K.; Kim, E.J.; Balázsi, G.; Rosner, M.R. Network of mutually repressive metastasis regulators can promote cell heterogeneity and metastatic transitions. Proc. Natl. Acad. Sci. USA 2014, 111, E364–E373. [Google Scholar] [CrossRef] [Green Version]
  51. Baritaki, S.; Huerta-Yepez, S.; Sahakyan, A.; Karagiannides, I.; Bakirtzi, K.; Jazirehi, A.; Bonavida, B. Mechanisms of nitric oxide-mediated inhibition of EMT in cancer. Cell Cycle 2010, 9, 4931–4940. [Google Scholar] [CrossRef] [Green Version]
  52. Yun, J.; Frankenberger, C.A.; Kuo, W.L.; Boelens, M.C.; Eves, E.M.; Cheng, N.; Liang, H.; Li, W.H.; Ishwaran, H.; Minn, A.J.; et al. Signalling pathway for RKIP and Let-7 regulates and predicts metastatic breast cancer. EMBO J. 2011, 30, 4500–4514. [Google Scholar] [CrossRef]
  53. Du, Y.; Zhu, H.C.; Liu, X.H.; Wang, L.; Ning, J.Z.; Xiao, C.C. MiR-543 Promotes Proliferation and Epithelial-Mesenchymal Transition in Prostate Cancer via Targeting RKIP. Cell. Physiol. Biochem. 2017, 41, 1135–1146. [Google Scholar] [CrossRef]
  54. de Castro, J.; Odeh, H.N.; Figy, C.; Yeung, M.L.; Trumbly, R.; Yeung, K.C. Regulation of RKIP expression in breast cancer cells by miRNAs. In Prognostic and Therapeutic Applications of RKIP in Cancer; Elsevier: Amsterdam, The Netherlands, 2020; pp. 139–146. [Google Scholar] [CrossRef]
  55. Du, Y.; Weng, X.D.; Wang, L.; Liu, X.H.; Zhu, H.C.; Guo, J.; Ning, J.Z.; Xiao, C.C. LncRNA XIST acts as a tumor suppressor in prostate cancer through sponging miR-23a to modulate RKIP expression. Oncotarget 2017, 8, 94358–94370. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  56. Fu, Z.; Kitagawa, Y.; Shen, R.; Shah, R.; Mehra, R.; Rhodes, D.; Keller, P.J.; Mizokami, A.; Dunn, R.; Chinnaiyan, A.M.; et al. Metastasis suppressor gene Raf kinase inhibitor protein (RKIP) is a novel prognostic marker in prostate cancer. Prostate 2006, 66, 248–256. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  57. Al-Mulla, F.; Bitar, M.S.; Al-Maghrebi, M.; Behbehani, A.I.; Al-Ali, W.; Rath, O.; Doyle, B.; Tan, K.Y.; Pitt, A.; Kolch, W. Raf Kinase Inhibitor Protein RKIP Enhances Signaling by Glycogen Synthase Kinase-3β. Cancer Res. 2011, 71, 1334–1343. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  58. Guo, W.; Dong, Z.; Lin, X.; Zhang, M.; Kuang, G.; Zhu, T. Decreased Expression and Aberrant Methylation of Raf Kinase Inhibitory Protein Gene in Esophageal Squamous Cell Carcinoma. Cancer Investig. 2012, 30, 703–711. [Google Scholar] [CrossRef]
  59. Kim, G.E.; Kim, N.I.; Lee, J.S.; Park, M.H.; Yoon, J.H. Reduced RKIP Expression is Associated With Breast Neoplastic Progression and is Correlated With Poor Outcomes and Aberrant Methylation in Breast Carcinoma. Appl. Immunohistochem. Mol. Morphol. 2017, 25, 467–474. [Google Scholar] [CrossRef]
  60. Minoo, P.; Baker, K.; Goswami, R.; Chong, G.; Foulkes, W.D.; Ruszkiewicz, A.R.; Barker, M.; Buchanan, D.; Young, J.; Jass, J.R. Extensive DNA methylation in normal colorectal mucosa in hyperplastic polyposis. Gut 2006, 55, 1467–1474. [Google Scholar] [CrossRef] [Green Version]
  61. Bonavida, B. Historical Perspectives of the Role of NO/NO Donors in Anti-tumor Activities: Acknowledging Dr. Keefer’s Pioneering Research. Crit. Rev. Oncog. 2021. [Google Scholar] [CrossRef]
  62. Fukumura, D.; Kashiwagi, S.; Jain, R.K. The role of nitric oxide in tumour progression. Nat. Rev. Cancer 2006, 6, 521–534. [Google Scholar] [CrossRef]
  63. Huerta-Yepez, S.; Vega, M.; Jazirehi, A.; Garban, H.; Hongo, F.; Cheng, G.; Bonavida, B. Nitric oxide sensitizes prostate carcinoma cell lines to TRAIL-mediated apoptosis via inactivation of NF-κB and inhibition of Bcl-xL expression. Oncogene 2004, 23, 4993–5003. [Google Scholar] [CrossRef] [Green Version]
  64. Bonavida, B. RKIP-Mediated Chemo-Immunosensitization of Resistant Cancer Cells via Disruption of the NF-κB/Snail/YY1/RKIP Resistance-Driver Loop. Crit. Rev. Oncog. 2014, 19, 431–445. [Google Scholar] [CrossRef]
  65. Rapozzi, V.; Pietra, E.D.; Zorzet, S.; Zacchigna, M.; Bonavida, B.; Xodo, L.E. Nitric oxide-mediated activity in anti-cancer photodynamic therapy. Nitric Oxide 2013, 30, 26–35. [Google Scholar] [CrossRef]
  66. Bonavida, B.; Baritaki, S. The Novel Role of Yin Yang 1 in the Regulation of Epithelial to Mesenchymal Transition in Cancer Via the Dysregulated NF-κB/Snail/YY1/RKIP/PTEN Circuitry. Crit. Rev. Oncog. 2011, 16, 211–226. [Google Scholar] [CrossRef] [PubMed]
  67. Pervin, S.; Singh, R.; Hernandez, E.; Wu, G.; Chaudhuri, G. Nitric Oxide in Physiologic Concentrations Targets the Translational Machinery to Increase the Proliferation of Human Breast Cancer Cells: Involvement of Mammalian Target of Rapamycin/eIF4E Pathway. Cancer Res. 2007, 67, 289–299. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  68. Pervin, S.; Tran, L.; Urman, R.; Braga, M.; Parveen, M.; Li, S.A.; Chaudhuri, G.; Singh, R. Oxidative stress specifically downregulates survivin to promote breast tumour formation. Br. J. Cancer 2013, 108, 848–858. [Google Scholar] [CrossRef] [Green Version]
  69. Rapozzi, V.; Pietra, E.D.; Bonavida, B. Dual roles of nitric oxide in the regulation of tumor cell response and resistance to photodynamic therapy. Redox Biol. 2015, 6, 311–317. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  70. Wink, D.A.; Ridnour, L.A.; Hussain, S.P.; Harris, C.C. The reemergence of nitric oxide and cancer. Nitric Oxide 2008, 19, 65–67. [Google Scholar] [CrossRef] [Green Version]
  71. Pasqualetti, G.; Ricciardi, S.; Mey, V.; Tacca, M.D.; Danesi, R. Synergistic cytotoxicity, inhibition of signal transduction pathways and pharmacogenetics of sorafenib and gemcitabine in human NSCLC cell lines. Lung Cancer 2011, 74, 197–205. [Google Scholar] [CrossRef]
  72. Giovannetti, E.; Labots, M.; Dekker, H.; Galvani, E.; Lind, J.; Sciarrillo, R.; Honeywell, R.; Smit, E.; Verheul, H.; Peters, G. Molecular Mechanisms and Modulation of Key Pathways Underlying the Synergistic Interaction of Sorafenib with Erlotinib in Non-Small-Cell-Lung Cancer (NSCLC) Cells. Curr. Pharm. Des. 2012, 19, 927–939. [Google Scholar] [CrossRef]
  73. Ricciardi, S.; Mey, V.; Nannizzi, S.; Pasqualetti, G.; Crea, F.; Tacca, M.D.; Danesi, R. Synergistic Cytotoxicity and Molecular Interaction on Drug Targets of Sorafenib and Gemcitabine in Human Pancreas Cancer Cells. Chemotherapy 2010, 56, 303–312. [Google Scholar] [CrossRef]
  74. Chatterjee, D.; Bai, Y.; Wang, Z.; Beach, S.; Mott, S.; Roy, R.; Braastad, C.; Sun, Y.; Mukhopadhyay, A.; Aggarwal, B.B.; et al. RKIP Sensitizes Prostate and Breast Cancer Cells to Drug-induced Apoptosis. J. Biol. Chem. 2004, 279, 17515–17523. [Google Scholar] [CrossRef] [Green Version]
  75. Yousuf, S.; Duan, M.; Moen, E.L.; Cross-Knorr, S.; Brilliant, K.; Bonavida, B.; LaValle, T.; Yeung, K.C.; Al-Mulla, F.; Chin, E.; et al. Raf Kinase Inhibitor Protein (RKIP) Blocks Signal Transducer and Activator of Transcription 3 (STAT3) Activation in Breast and Prostate Cancer. PLoS ONE 2014, 9, e92478. [Google Scholar] [CrossRef]
  76. Hampl, V.; Tristani-Firouzi, M.; Hutsell, T.C.; Archer, S.L. Nebulized nitric oxide/nucleophile adduct reduces chronic pulmonary hypertension. Cardiovasc. Res. 1996, 31, 55–62. [Google Scholar] [CrossRef]
  77. Wishart, D.S.; Feunang, Y.D.; Guo, A.C.; Lo, E.J.; Marcu, A.; Grant, J.R.; Sajed, T.; Johnson, D.; Li, C.; Sayeeda, Z.; et al. DrugBank 5.0: A major update to the DrugBank database for 2018. Nucleic Acids Res. 2017, 46, D1074–D1082. [Google Scholar] [CrossRef]
  78. Walker, E.J.; Rosenberg, S.A.; Wands, J.R.; Kim, M. Role of Raf Kinase Inhibitor Protein in Hepatocellular Carcinoma. For. Immunopathol. Dis. Therap. 2011, 2, 195–204. [Google Scholar] [CrossRef] [PubMed]
  79. Wu, X.; Yang, Y.; Xu, Z.; Li, J.; Yang, B.; Feng, N.; Zhang, Y.; Wang, S. Raf kinase inhibitor protein mediated signaling inhibits invasion and metastasis of hepatocellular carcinoma. Biochim. Biophys. Acta 2016, 1860, 384–391. [Google Scholar] [CrossRef] [PubMed]
  80. Lei, X.; Chang, L.; Ye, W.; Jiang, C.; Zhang, Z. Raf kinase inhibitor protein (RKIP) inhibits the cell migration and invasion in human glioma cell lines in vitro. Int. J. Clin. Exp. Pathol. 2015, 8, 14214–14220. [Google Scholar] [PubMed]
  81. Zhao, D.; Ma, J.; Shi, J.; Cheng, L.; Li, F.; Jiang, X.; Jiang, H. Raf kinase inhibitor protein inhibits esophageal cancer cell invasion through downregulation of matrix metalloproteinase expression. Oncol. Rep. 2013, 30, 304–312. [Google Scholar] [CrossRef]
  82. Hao, C.; Wei, S.; Tong, Z.; Li, S.; Shi, Y.; Wang, X.; hua Zhu, Z. The effects of RKIP gene expression on the biological characteristics of human triple-negative breast cancer cells in vitro. Tumour Biol. 2012, 33, 1159–1167. [Google Scholar] [CrossRef]
  83. Zou, Q.; Wu, H.; Fu, F.; Yi, W.; Pei, L.; Zhou, M. RKIP suppresses the proliferation and metastasis of breast cancer cell lines through up-regulation of miR-185 targeting HMGA2. Arch. Biochem. Biophys. 2016, 610, 25–32. [Google Scholar] [CrossRef] [PubMed]
  84. Lin, K.; Baritaki, S.; Militello, L.; Malaponte, G.; Bevelacqua, Y.; Bonavida, B. The Role of B-RAF Mutations in Melanoma and the Induction of EMT via Dysregulation of the NF- B/Snail/RKIP/PTEN Circuit. Genes Cancer 2010, 1, 409–420. [Google Scholar] [CrossRef] [PubMed]
  85. Schuierer, M.M.; Bataille, F.; Hagan, S.; Kolch, W.; Bosserhoff, A.K. Reduction in Raf Kinase Inhibitor Protein Expression Is Associated with Increased Ras-Extracellular Signal-Regulated Kinase Signaling in Melanoma Cell Lines. Cancer Res. 2004, 64, 5186–5192. [Google Scholar] [CrossRef] [Green Version]
  86. Minoo, P.; Zlobec, I.; Baker, K.; Tornillo, L.; Terracciano, L.; Jass, J.R.; Lugli, A. Loss of Raf-1 Kinase Inhibitor Protein Expression Is Associated With Tumor Progression and Metastasis in Colorectal Cancer. Am. J. Clin. Pathol. 2007, 127, 820–827. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  87. Lee, T.Y.; Liu, C.L.; Chang, Y.C.; Nieh, S.; Lin, Y.S.; Jao, S.W.; Chen, S.F.; Liu, T.Y. Increased chemoresistance via Snail-Raf kinase inhibitor protein signaling in colorectal cancer in response to a nicotine derivative. Oncotarget 2016, 7, 23512–23520. [Google Scholar] [CrossRef]
  88. Morozova, N.; Zinovyev, A.; Nonne, N.; Pritchard, L.L.; Gorban, A.N.; Harel-Bellan, A. Kinetic signatures of microRNA modes of action. RNA 2012, 18, 1635–1655. [Google Scholar] [CrossRef] [Green Version]
  89. Hou, W.; Shan, Y.; Zheng, J.; Lambrecht, R.W.; Donohue, S.E.; Bonkovsky, H.L. Zinc mesoporphyrin induces rapid and marked degradation of the transcription factor bach1 and up-regulates HO-1. Biochim. Biophys. Acta Gene Regul. Mech. 2008, 1779, 195–203. [Google Scholar] [CrossRef] [Green Version]
  90. Shahrezaei, V.; Swain, P.S. Analytical distributions for stochastic gene expression. Proc. Natl. Acad. Sci. USA 2008, 105, 17256–17261. [Google Scholar] [CrossRef] [Green Version]
  91. Yan, X.; Ehnert, S.; Culmes, M.; Bachmann, A.; Seeliger, C.; Schyschka, L.; Wang, Z.; Rahmanian-Schwarz, A.; Stöckle, U.; Sousa, P.A.D.; et al. 5-Azacytidine Improves the Osteogenic Differentiation Potential of Aged Human Adipose-Derived Mesenchymal Stem Cells by DNA Demethylation. PLoS ONE 2014, 9, e90846. [Google Scholar] [CrossRef] [Green Version]
  92. Deng, X.; Liu, Z.; Li, X.; Zhou, Y.; Hu, Z. Generation of new hair cells by DNA methyltransferase (Dnmt) inhibitor 5-azacytidine in a chemically-deafened mouse model. Sci. Rep. 2019, 9, 7997. [Google Scholar] [CrossRef] [PubMed]
  93. Hays, E.; Bonavida, B. Nitric Oxide-Mediated Enhancement and Reversal of Resistance of Anticancer Therapies. Antioxidants 2019, 8, 407. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  94. Arfken, G.B.; Weber, H.J. Mathematical Methods for Physicists, 6th ed.; Elsevier: Amsterdam, The Netherlands, 2005. [Google Scholar]
  95. Innocentini, G.C.P.; Hornos, J.E.M. Modeling stochastic gene expression under repression. J. Math. Biol. 2007, 55, 413–431. [Google Scholar] [CrossRef] [Green Version]
  96. Abramowitz, M.; Stegun, I.A. (Eds.) Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables, 10th ed.; National Bureau of Standards Applied Mathematics Series; U.S. Department of Commerce, National Bureau of Standards: Washington, DC, USA, 1972; Volume 55.
Figure 1. Drugs aiming at 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 on 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 a single or combination drugs following a dose agenda.
Figure 1. Drugs aiming at 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 on 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 a single or combination drugs following a dose agenda.
Cancers 14 00633 g001
Figure 2. Probability distributions. The red lines indicate pre-treatment probabilities, while the green lines indicate the aimed probabilities governing the post-treatment number of gene products. The post-treatment distribution parameters were set aiming to increase the average number of transcripts to be ≈100. We disregarded the effect of the new parameter values on the variances, i.e., the heterogeneity of treatment response. We consider five initial conditions for the gene switching speed, ϵ 0 = ( 0.1, 1 , 2 , 10 , 10 ) , arranged in rows and indicated by (A1A3), (B1B3), (C1C3), (D1D3) and (E1E3) in the labels of respective graphs. (D,E) have the same value for ϵ 0 but represent different because of the differing values of A 0 . The distributions for the treatment designs are arranged in columns: (1) DETANONOate aiming the f rate, (2) 5-AzaC aiming the k rate, and (3) both drugs aiming both f, and k rates, simultaneously.
Figure 2. Probability distributions. The red lines indicate pre-treatment probabilities, while the green lines indicate the aimed probabilities governing the post-treatment number of gene products. The post-treatment distribution parameters were set aiming to increase the average number of transcripts to be ≈100. We disregarded the effect of the new parameter values on the variances, i.e., the heterogeneity of treatment response. We consider five initial conditions for the gene switching speed, ϵ 0 = ( 0.1, 1 , 2 , 10 , 10 ) , arranged in rows and indicated by (A1A3), (B1B3), (C1C3), (D1D3) and (E1E3) in the labels of respective graphs. (D,E) have the same value for ϵ 0 but represent different because of the differing values of A 0 . The distributions for the treatment designs are arranged in columns: (1) DETANONOate aiming the f rate, (2) 5-AzaC aiming the k rate, and (3) both drugs aiming both f, and k rates, simultaneously.
Cancers 14 00633 g002
Figure 3. The dynamics of the number of RKIP mRNA (green solid line) and their standard deviation (red solid line) along the time for treatment aiming f kinetic rate. The black dashed line at 80 indicates an arbitrary threshold on the number of RKIP mRNA. (A1E1) show single doses of drug with different progressive values to ϵ 0 (shown within the respective graph). (F1J1) show the same initial conditions that (A1E1) for 5 doses with time interval of 10 h (indicated by arrows). All doses are in maximum tolerance ( ξ a = 1 ).
Figure 3. The dynamics of the number of RKIP mRNA (green solid line) and their standard deviation (red solid line) along the time for treatment aiming f kinetic rate. The black dashed line at 80 indicates an arbitrary threshold on the number of RKIP mRNA. (A1E1) show single doses of drug with different progressive values to ϵ 0 (shown within the respective graph). (F1J1) show the same initial conditions that (A1E1) for 5 doses with time interval of 10 h (indicated by arrows). All doses are in maximum tolerance ( ξ a = 1 ).
Cancers 14 00633 g003
Figure 4. The dynamics of the number of RKIP mRNA (green solid line) and their standard deviation (red solid line) along the time for treatment aiming k kinetic rate. The black dashed line at 80 indicates an arbitrary threshold on the number of RKIP mRNA. (A2E2) show single doses of drugs with different progressive values to ϵ 0 (shown within the respective graph). (F2J2) show the same initial conditions that (A2E2) for 10 doses with time interval of 4 h (indicated by arrows). All doses are in maximum tolerance ( ξ b = 1 ).
Figure 4. The dynamics of the number of RKIP mRNA (green solid line) and their standard deviation (red solid line) along the time for treatment aiming k kinetic rate. The black dashed line at 80 indicates an arbitrary threshold on the number of RKIP mRNA. (A2E2) show single doses of drugs with different progressive values to ϵ 0 (shown within the respective graph). (F2J2) show the same initial conditions that (A2E2) for 10 doses with time interval of 4 h (indicated by arrows). All doses are in maximum tolerance ( ξ b = 1 ).
Cancers 14 00633 g004
Figure 5. The dynamics of the number of RKIP mRNA along the time for treatment with two drugs that aim at f and k kinetic rates are shown in green solid lines, and their standard deviations are shown in red solid lines. The threshold for RKIP is 80, and it is shown within graphs as a black dashed line. (A3E3) show the single doses of both drugs for different progressive values to ϵ 0 (indicated within each graph). (F3J3) show the same initial conditions that (A3E3) for multiple fractional doses of drug that target f (or k) with time interval of 10 h (or 4 h), in order, with 8 (or 20) doses in (F3) and 6 (or 15) doses in (G3J3). Arrows indicate the agenda, with time moments and fractions ( ξ a and ξ b ) for each drug, DETANONOate (aiming f) in blue and 5-AzaC (aiming k) in magenta. The numbers on the larger arrows indicate the fraction of the dose applied at that time and those that follow (arrows without numbers). The agenda is also shown in Table 1.
Figure 5. The dynamics of the number of RKIP mRNA along the time for treatment with two drugs that aim at f and k kinetic rates are shown in green solid lines, and their standard deviations are shown in red solid lines. The threshold for RKIP is 80, and it is shown within graphs as a black dashed line. (A3E3) show the single doses of both drugs for different progressive values to ϵ 0 (indicated within each graph). (F3J3) show the same initial conditions that (A3E3) for multiple fractional doses of drug that target f (or k) with time interval of 10 h (or 4 h), in order, with 8 (or 20) doses in (F3) and 6 (or 15) doses in (G3J3). Arrows indicate the agenda, with time moments and fractions ( ξ a and ξ b ) for each drug, DETANONOate (aiming f) in blue and 5-AzaC (aiming k) in magenta. The numbers on the larger arrows indicate the fraction of the dose applied at that time and those that follow (arrows without numbers). The agenda is also shown in Table 1.
Cancers 14 00633 g005
Figure 6. Effectiveness of treatment agendas measured by the average values n σ ¯ (top row) and 2 σ ¯ (bottom row) during the interval from 60 to 80 h after application of the first dose is presented. The vertical bars on the right give color code denoting the average values. In each square of a heatmap, we indicate the value of n σ ¯ (top row) and 2 σ ¯ (bottom row) as a result of the fractional dose effect ξ a and ξ b on rates f and k, respectively. The analysis was performed for the five pre-treatment conditions identified by the switching speed value ϵ 0 labeling each column. For each pair ( ξ a , ξ b ), we simulated the application of DETANONOate (and 5-AzaC) every 10 h (and 4 h). For a treatment aiming RKIP gene, dose reduction to minimize cytotoxic effects for each ϵ 0 consists on keeping ( ξ a , ξ b ) within regions with n σ ¯ above Graphs (AE) and 2 σ ¯ below Graphs (FJ) specific thresholds. The yellow regions of the heatmaps indicate the more desirable domains evaluated in terms of n σ ¯ and 2 σ ¯ . The sequence of red arrows in each graph corresponds to the effect of dose fraction combinations in the agendas in Figure 5 and Table 1, for each pre-treatment condition.
Figure 6. Effectiveness of treatment agendas measured by the average values n σ ¯ (top row) and 2 σ ¯ (bottom row) during the interval from 60 to 80 h after application of the first dose is presented. The vertical bars on the right give color code denoting the average values. In each square of a heatmap, we indicate the value of n σ ¯ (top row) and 2 σ ¯ (bottom row) as a result of the fractional dose effect ξ a and ξ b on rates f and k, respectively. The analysis was performed for the five pre-treatment conditions identified by the switching speed value ϵ 0 labeling each column. For each pair ( ξ a , ξ b ), we simulated the application of DETANONOate (and 5-AzaC) every 10 h (and 4 h). For a treatment aiming RKIP gene, dose reduction to minimize cytotoxic effects for each ϵ 0 consists on keeping ( ξ a , ξ b ) within regions with n σ ¯ above Graphs (AE) and 2 σ ¯ below Graphs (FJ) specific thresholds. The yellow regions of the heatmaps indicate the more desirable domains evaluated in terms of n σ ¯ and 2 σ ¯ . The sequence of red arrows in each graph corresponds to the effect of dose fraction combinations in the agendas in Figure 5 and Table 1, for each pre-treatment condition.
Cancers 14 00633 g006
Figure 7. The enhanced treatment designs resulting from distributions of Figure 2A2,B3,C2,D3,E3. We simulate a hypothetical cocktail of drugs changing all rates ( f , k , h , ρ ) . The first row of graphs shows the pre- (red lines) and post-(green lines) treatment probability distributions governing the mRNA numbers. An asterisk is added to the labels to indicate that we are simulating the enhanced treatment design. The average number n for all pre-treatment distributions is 10 and for each post-treatment distribution is ( 103 , 88 , 63 , 63 , 44 ) following graphs from left to right, respectively. The second row shows the respective dynamic of responses to enhanced treatment design in (A4E4). Green lines indicate the average numbers of RKIP mRNAs, n . Red lines show the one standard deviation around the average, n ± σ . Black dashed lines at n = 80 indicate the threshold separating the mRNA expression levels in a healthy (above line) and cancer (below) cell. The parameters of post-treatment distributions and the agenda of doses of each enhanced treatment design are shown in Appendix B.
Figure 7. The enhanced treatment designs resulting from distributions of Figure 2A2,B3,C2,D3,E3. We simulate a hypothetical cocktail of drugs changing all rates ( f , k , h , ρ ) . The first row of graphs shows the pre- (red lines) and post-(green lines) treatment probability distributions governing the mRNA numbers. An asterisk is added to the labels to indicate that we are simulating the enhanced treatment design. The average number n for all pre-treatment distributions is 10 and for each post-treatment distribution is ( 103 , 88 , 63 , 63 , 44 ) following graphs from left to right, respectively. The second row shows the respective dynamic of responses to enhanced treatment design in (A4E4). Green lines indicate the average numbers of RKIP mRNAs, n . Red lines show the one standard deviation around the average, n ± σ . Black dashed lines at n = 80 indicate the threshold separating the mRNA expression levels in a healthy (above line) and cancer (below) cell. The parameters of post-treatment distributions and the agenda of doses of each enhanced treatment design are shown in Appendix B.
Cancers 14 00633 g007
Table 1. The fractions ξ a and ξ b of the maximal tolerated dose for treatment agendas of the pre-treatment conditions in Figure 5.
Table 1. The fractions ξ a and ξ b of the maximal tolerated dose for treatment agendas of the pre-treatment conditions in Figure 5.
GraphSequence of Number of Doses × Fraction forCumulative Reduction in
ξ a ξ b ξ a ξ b
F3 3 × 0.9 ; 4 × 0.8 ; 1 × 0.5 5 × 0.9 ; 10 × 0.8 ; 5 × 0.7 20 % 20 %
G3 3 × 0.8 ; 1 × 0.7 ; 2 × 0.5 5 × 0.8 ; 10 × 0.75 32 % 23 %
H3 3 × 0.8 ; 1 × 0.7 ; 2 × 0.5 5 × 0.8 ; 10 × 0.75 32 % 23 %
I3 4 × 0.7 ; 1 × 0.6 ; 1 × 0.5 15 × 0.7 35 % 30 %
J3 3 × 0.7 ; 1 × 0.65 ; 2 × 0.6 5 × 0.7 ; 5 × 0.65 ; 5 × 0.6 29 % 35 %
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Giovanini, G.; Barros, L.R.C.; Gama, L.R.; Tortelli, T.C., Jr.; Ramos, A.F. A Stochastic Binary Model for the Regulation of Gene Expression to Investigate Responses to Gene Therapy. Cancers 2022, 14, 633. https://doi.org/10.3390/cancers14030633

AMA Style

Giovanini G, Barros LRC, Gama LR, Tortelli TC Jr., Ramos AF. A Stochastic Binary Model for the Regulation of Gene Expression to Investigate Responses to Gene Therapy. Cancers. 2022; 14(3):633. https://doi.org/10.3390/cancers14030633

Chicago/Turabian Style

Giovanini, Guilherme, Luciana R. C. Barros, Leonardo R. Gama, Tharcisio C. Tortelli, Jr., and Alexandre F. Ramos. 2022. "A Stochastic Binary Model for the Regulation of Gene Expression to Investigate Responses to Gene Therapy" Cancers 14, no. 3: 633. https://doi.org/10.3390/cancers14030633

APA Style

Giovanini, G., Barros, L. R. C., Gama, L. R., Tortelli, T. C., Jr., & Ramos, A. F. (2022). A Stochastic Binary Model for the Regulation of Gene Expression to Investigate Responses to Gene Therapy. Cancers, 14(3), 633. https://doi.org/10.3390/cancers14030633

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop