Stochastic modelling and simulation of PTEN regulatory networks with miRNAs and ceRNAs

In this work, three genetic regulatory networks are considered, that model the post–transcriptional regulation of the PTEN onco–suppressor gene, mediated by microRNAs and competitive endogenous RNAs, in glioblastoma multiforme, the most severe of brain tumours. We simulate solutions of the resulting stochastic differential systems and discuss the effects of this miRNA–fashioned regulation on PTEN expression.


Introduction
In the human body, cells contain the approximately twenty thousand genes encoded in the human genome. Genes get selectively expressed (or turned on) within each cell, thus determining its type, shape and behaviour. A gene regulatory network (GRN) is a network inferred from gene expression data [11], so it is a tool for describing gene-gene (as well as potential protein-protein) interactions and the way genes can turn each other on and off. In other words, a GRN provides information about the regulatory interactions between regulators and their potential targets, interactions that can be incredibly complicated, even in a network of very few genes, allowing a cell to respond to signals in various complex ways [11].
GRN modelling represents, therefore, an important postgenomic-era challenge, and it has indeed gained more and more interest, also in computational bio-mathematics. Several methods have been proposed for estimating gene networks from gene expression data, where both stochastic and deterministic approaches are pursued [4,5,22,28].
In this work, the focus is on the regulation of PTEN (Phosphatase and TENsin homolog) gene expression specifically in the brain, considering transcriptional and post-transcriptional interactions, and excluding translation in the PTEN protein, for which we refer to [8]. PTEN is found at locus 10q23.31 on chromosome 10, and includes instructions to produce a protein that helps govern the cell life-death cycle by preventing cells from growing and multiplying too rapidly and uncontrollably. PTEN is thus an onco-suppressor, but it is also one of the most frequently mutated or silenced in human cancer [9]. Loss of chromosomal segment 10q is indeed proved to be linked to glioblastoma multiforme (GBM), a very aggressive brain tumour [21,29]. Even in the absence of genetic loss or mutation, therefore in a situation where PTEN levels should be normal, PTEN is found to be low in about 70 % of GBM cases, which prompted research on PTEN regulation in an attempt to justify this unexplained low expression.
In the GRNs we consider, the crucial factors involved are short non-coding RNAs (microRNAs or miRNAs) and competitive endogenous RNAs (ceRNAs) [25]. To understand their role in gene expression control, it is useful to synthesise it as a process consisting of a sequence of steps, involving transcription of DNA to mRNA, namely an RNA intermediary or messenger RNA, and translation of mRNA into the final functional product, the protein. After transcription, miRNAs negatively regulate expression of their target mRNAs, where one miRNA can regulate many mRNAs, while the same mRNA can be regulated by a number of miRNAs. As a consequence of this down-regulation, a miRNA can be an onco-suppressor or an onco-promoter. The activity of miRNAs can be influenced by the presence of ceRNAs: a miRNA controls expression of a gene by binding to its mRNA, and ceRNAs compete for miRNA binding sites. The result is a complex network, in which miRNAs and ceRNAs interact and affect the mRNA level of the gene considered (PTEN in our case), available to be translated into the final protein.
The importance of PTEN ceRNA-mediated regulation in tumour cells is studied in [17,25,26], but a clear quantitative understanding of its miRNA-mediated regulation is still lacking [12]. One attempt is represented by the mass-action mathematical deterministic model in [1], used to determine optimal conditions for ceRNAs activity in silico. A stochastic and discrete model can be found in [8], where translation into PTEN protein and transcriptional and translational delays are also considered. A complex network of a miRNA pool acting in the brain, including corresponding ceRNAs competing with PTEN for the same miRNAs, is introduced and analysed numerically in [7], where the continuous and deterministic setting, modelled by ordinary differential equations, was used to overcome computational and time complexity.
In the present article, we model the effects of three of the miRNAs included in the pool studied in [7], one at a time, so to make clearer their specific contribution within the network [20]. Given the reduced computational complexity of the reacting systems, the more realistic stochastic and continuous (rather than deterministic) modelling regime is employable here.

GRNs modelling and simulation regimes
When simulating a Gene Regulatory Newtwork (GRN), three modelling regimes can be considered: stochastic and discrete, stochastic and continuous, deterministic and continuous; additional complexity exists, related to simultaneously considering temporal and spatial effects, which is not dealt with in this article [19].
The discrete stochastic regime (also called slow) is characterised by the use of a Stochastic Simulation Algorithm (SSA), the first of which is due to Gillespie [13]. It is an essentially exact procedure for numerically simulating the evolution of a set of chemical reactions in a well-stirred and homogeneous chemical reaction system, taking into account inherent randomness.
If X 1 , . . . , X n s represent molecular species, interacting through R 1 , . . . , R n r reaction channels, then the system evolution is described by a discrete non-linear Markov process in which a vector X(t), holding numbers (integer values) of the molecular species at time t, evolves through time.
The state vector X(t) is thus a discrete jump Markov process, whose time evolution equation, the so-called Chemical Master Equation (CME) reminded in (1) below, describes the probability P(X, t | X 0 , t 0 ) that, given an initial condition X(t 0 ) = X 0 , then X(t) = X . Now, any set of n r chemical reactions is uniquely characterised by stoichiometric vectors ν 1 , . . . , ν n r and by propensity functions a 1 (X(t)), . . . , a n r (X(t)). Vectors ν j , each of dimension n s , represent the update of the numbers of molecules in the system if R 1 , . . . , R n r respectively occur, while functions a j = a j (X(t)) describe the relative probabilities of R 1 , . . . , R n r occurring, again respectively.
Hence, the CME takes the following form: In general, (1) is a discrete parabolic partial differential equation too difficult to solve, either analytically or numerically: other techniques are therefore needed to simulate X(t) .
SSAs, which are rigorously based on the same microphysical premise that underlies the CME, can be employed when there are low molecular numbers and intrinsic, or internal, noise cannot be ignored; the latter is linked to the fact that it is unknown which reaction is firing, and at what time it is firing.
An SSA yields a more realistic representation of a system evolution, compared to the deterministic Reaction Rate Equation (RRE), i.e. system (2) of ordinary differential equations that ignores the presence of noise: where a = (a 1 , . . . , a n r ) T is the propensity vector, while S = ν k j is the stoichiometric matrix of dimensions n s × n r , whose columns are the stoichiometric vectors. The RRE intervenes when the chemical reaction system is too large to be treated within the discrete stochastic setting. The SSA, in fact, remains a computationally demanding approach, though many authors have attempted to enhance its performance [6,16]; this limits its applicability, in particular to the large reaction networks required for modelling realistic gene networks.
In such a situation, as said, the RRE can come into play, namely solution is sought within the third regime, deterministic and continuous, aimed at describing, in an averaged manner, the behaviour of a system in which the number of molecules is so large that we can speak of concentrations of reactants and products. This averaged behaviour can be modelled as an Initial Value Problem formed by (2) together with an initial condition The third regime (also named fast) is, indeed, the standard one in chemical kinetics where the Law of Mass Action applies, i.e. the solution to the RRE satisfies this law, and is mainly used to model complex GRNs [7], when the first and second regimes fail computationally, or when the deterministic component dominates the noise terms. Note, though, that the RRE is not appropriate if the molecular population of some critical reactant species is so small that microscopic fluctuations can produce macroscopic effects, which is essentially true for genetic/enzymatic reactions in living cells [2,10,15].
In this article, we consider the second regime, where noise effects are still important, but continuity arguments hold up; this continuous stochastic regime is also referred to as intermediate regime, since it represents a good compromise between accuracy of SSAs and fastness of RRE solvers.
By applying the Central Limit Theorem and matching the first two moments of the CME, it is possible to write a stochastic ordinary differential equation (SODE) that describes the evolution of X(t) [14]. This equation is the so-called Chemical Langevin Equation (CLE) and takes the form: where matrix B = S · a · S T 1/2 and the diagonal matrix a = diag(a) have dimensions n s × n s and n r × n r , respectively, while vector W = W(t) = is formed by independent standard Wiener processes, each with zero drift and unit volatility.
Equation (3) is an Itô SODE that models the continuous stochastic evolution equation and takes into account intrinsic noise [14].
General classes of methods can be used to solve (3) in non-stiff cases. The simplest is the Euler-Maruyama method, where the stochastic increments, at each time step, take the form: These Wiener increments are distributed according to a normal random variable of zero mean and variance given by their length h . The Euler-Maruyama method converges with strong order 1/2 to the Itô solution of (3) [3,18,27].
Before leaving this section, we recall that both SSA and (3) admit solutions converging, in the limit of large numbers of reactants, to the solution of (2), thus satisfying the Law of Mass Action [14].

miRNAs
Three particular miRNAs, involved in the onset of brain cancer, are considered in Sects. 3.1, 3.2, and 3.3. For each of them, the interaction is modelled with the PTEN onco-suppressor gene, as well as other genes concurrent with PTEN, namely the ceRNAs. In particular, for miRNA 214 , the PTENP1 pseudo-gene acts as ceRNA; for miRNA 19b , CNOT6L gene and PTENP1 pseudo-gene act as ceRNAs; finally, for miRNA 17 , VAPA and CNOT6L genes and PTENP1 pseudo-gene act as ceRNAs.

miRNA 214
We consider a set of n r chemical reactions, monomolecular or bimolecular with distinct reactants, where n s is the total number of X-species that interact through such reactions. Set (4) represents a realistic example of a system of ten reactions involving seven X-species.  By introducing differential variables X k = X k (t), each representing the k-th Xspecies observed and its biological state at time t, the reaction set can be translated into a system of n s differential equations in n s variables. Table 1 shows the correspondence between X-species and differential variables, so that (4) becomes (5). We seek its solution X(t) = (X 1 (t), . . . , X n s (t)) T , namely the state vector whose evolution in time t is that of the differential system under consideration. To ease the notation, in the following we often set X = X(t).
For each j-th reaction, a propensity function a j = a j (X(t)) is defined, to express the link between the reactants involved and the probability c j (per time unit) that a randomly chosen combination of these reactants will indeed react. In the numerical simulations presented in Sect. 4, the c j values are inferred according to what is known in the literature, for which we refer the reader to [7] and references therein, and in the GeneCard database [24].
The propensity functions corresponding to reactions (5) are shown in equalities (6), where dependence on time t has been temporarily discarded, again to simplify the notation. a j = c j X j j = 1 . . . 6 , a 7 = c 7 X 4 , a 8 = c 8 X 2 X 7 , a 9 = c 9 X 6 X 7 , a 10 = c 10 X 7 .
In addition to the propensity vector a = (a 1 , . . . , a n r ) T , the other information needed, to set up the Itô process related to a chemical reaction system, is its stoichiometry.
As seen in Sect. 2, in fact, propensity and stoichiometry values are the ingredients to build the SODE (3) associated with the set of reactions under observation: vector a and the n s × n r matrix S = ν k j determine, in fact, both deterministic and stochastic components in Eq. (3). Table 2 Stoichiometric matrix S = ν k j for reactions (5) ν 1 ν 2 ν 3 ν 4 ν 5 ν 6 ν 7 ν 8 ν 9 ν 10 Table 3 Initial conditions X k (t 0 ) = X k,0 , with t 0 = 0, for reactions (5) k X k,0 k X k,0 k X k,0 k X k,0  Now, for each j-th reaction, the corresponding column of S is: Hence, in the case of system (5), the S-matrix is as in Table 2. Tables 3 and 4 contain, respectively, initial conditions and c j coefficients employed, in the numerical simulation Sect. 4, to set and solve the SODE associated with reactions (5).

miRNA 19b
Set (8) represents a realistic example of a thirteen reaction system involving nine Xspecies. This time, the chemical reactions are written directly in terms of differential variables. Table 5 shows the correspondence between X-species and differential variables in (8).
Tables 6 and 7 contain, respectively, initial conditions and c j coefficients employed, in Sect. 4, to set and solve the SODE associated with reactions (8). Table 6 Initial conditions X k (t 0 ) = X k,0 , with t 0 = 0, for reactions (8)

miRNA 17
Set (11) represents a realistic example of a system of sixteen reactions involving eleven X-species, where the correspondence between X-species and differential variables is given in Table 8.
Tables 9 and 10 contain, respectively, initial conditions and c j coefficients employed, in Sect. 4, to set and solve the SODE associated with reactions (11).  Table 9 Initial conditions X k (t 0 ) = X k,0 , with t 0 = 0, for reactions (11) k X k,0 k X k,0 k X k,0 k X k,0

Numerical simulations
In this section, systems (5), (8) and (11)   All experiments were carried out on a DELL XPS notebook computer, with 7th generation Intel Core i7 processor at 2.8 GHz, and 16 GB RAM at 2400 MHz, running under Ubuntu 20.04.4 LTS operating system. The problem solving environment of Mathematica, version 13, is employed to set up each model, integrate the resulting SODE, and plot its solution [23]. Timing results, averaged over ten runs for each experiment, show that the computational time taken by the definition of the Itô process is around seven seconds, most of which is spent to form the matrix B in (3), while its integration requires about six seconds.
In all experiments, the focus is on the interaction between PTEN and miRNA, that is X 2 and X 4 . This is done to establish whether down-regulation of X 2 or X 4 implies over-expression of the other. Furthermore, the evolution is plotted along five temporal intervals, to cover short-and long-term behaviours. Figure 1 illustrates the solution for reactions (5). Simulations show that miRNA 214 remains expressed at a very low costant level, so it does not affect the PTEN regulation. In fact, after a short initial phase during which there is a fairly dramatic increase of Fig. 3 Plots of species X 2 and X 4 , respectively PTEN and miRNA 17 , for reactions (11) over twenty minutes (top), one day (centre left), four days (centre right), one month (bottom left), two months (bottom right). miRNA 17 rapidly grows and, within a very short time, induces a dramatic decay in the PTEN level PTEN, its value settles around a constant mean value. This is exactly what we expect in a physiological, non-pathological, situation within a cell. In particular, it can be observed that the mean level of miRNA 214 remains close to that of its initial condition of order 10 0 , while that of PTEN oscillates around a value at an order of magnitude 10 2 . Figure 2 illustrates the solution for reactions (8). Simulations with miRNA 19b show that its expression grows, initially and rapidly, up to a value of about 150, and then it stabilises, reaching a situation similar to the miRNA 214 case. In fact, even if in this case both PTEN and miRNA 19b oscillate around a mean value at the same order of magnitude 10 2 , miRNA 19b has, once again, no effect on PTEN expression. Figure 3 illustrates the solution for reactions (11). Simulations with miRNA 17 show how it can cause a situation that is no longer physiological inside a cell, but rather pathological. Here, in fact, miRNA 17 grows dramatically to a mean value of order of magnitude 10 3 , large enough to induce the same dramatic down-regulation of PTEN. Due to this highly reduced function of PTEN as an onco-suppressor, the cell is greatly exposed to tumour occurence.

Conclusions
This article investigates the impact of three different miRNAs on the regulation of the PTEN gene, an onco-suppressor whose low expression is observed in the majority of cases of glioblastoma multiforme, one of the most aggressive brain tumours. We considered miRNA 214 , miRNA 19b and miRNA 17 , each targeting different genes other than PTEN, called ceRNAs. The three resulting GRNs are simulated in a stochastic and continuous setting, via systems of stochastic ordinary differential equations derived from jump Markov processes, modelling appropriate chemical reactions. Our simulations show that, unless miRNAs grow in such a way to exceed PTEN levels by at least one order of magnitude, PTEN preserves its physiological oscillatory behaviour, around a mean value of order 10 2 . The above situation occurs with miRNA 214 and miRNA 19b . In the case of miRNA 17 , simulations show that it rapidly becomes over-expressed, while PTEN is dramatically down-regulated, leaving the cell much more unprotected from tumour genesis. A more faithful adherence to reality could be obtained by considering the stochastic and discrete processes underlying the GRNs models. This choice is currently under investigation, leading to the use of a stochastic simulation algorithm and incorporating delays in the reactions that represent transcription of DNA to RNA, both for PTEN and ceRNAs. appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.