Single-molecule modeling of mRNA degradation by miRNA: Lessons from data

Recent experimental results on the effect of miRNA on the decay of its target mRNA have been analyzed against a previously hypothesized single molecule degradation pathway. According to that hypothesis, the silencing complex (miRISC) first interacts with its target mRNA and then recruits the protein complexes associated with NOT1 and PAN3 to trigger deadenylation (and subsequent degradation) of the target mRNA. Our analysis of the experimental decay patterns allowed us to refine the structure of the degradation pathways at the single molecule level. Surprisingly, we found that if the previously hypothesized network was correct, only about 7% of the target mRNA would be regulated by the miRNA mechanism, which is inconsistent with the available knowledge. Based on systematic data analysis, we propose the alternative hypothesis that NOT1 interacts with miRISC before binding to the target mRNA. Moreover, we show that when miRISC binds alone to the target mRNA, the mRNA is degraded more slowly, probably through a deadenylation-independent pathway. The new biochemical pathway proposed here both fits the data and paves the way for new experimental work to identify new interactions.


From single-molecule pathways to decay patterns
We model mRNA decay as a single molecule stochastic process. The process starts from an initial state 0, which represents the mRNA in its initial condition, and terminates in the degradation state n. The initial and the degradation states are defined here based on the time point from which the mRNA is experimentally detected until it is not detected anymore, respectively, in the context of the experimental technique used in [1]. Each realization of the stochastic process, which starts from state 0 to the degradation state through the network of states represents the life of a single mRNA molecule. Therefore, the distribution of the random time T elapsing from state 0 to the degradation state n can be identified with the distribution of the lifetimes of the mRNA molecules. The state space σ = {0, . . . , n} is thus made of n transient states and one absorbing state n (a generalization to more absorbing states is straightforward, see [2]).
In our modeling framework each state transition describes the occurrence of what we define a first-order biochemical event. This definition includes three possible biochemical scenarios: (i) an elementary first-order chemical reaction (ii) a pseudo first-order biochemical reaction: the actual reaction order is > 1 but the concentration all of the reactants except for one are at the steady state level. As a result, this reaction behaves as a "true" first-order reaction (iii) an apparent first-order reaction: a complex chain of reactions with a single rate limiting step, thus exhibiting first-order behavior In the context of the mRNA degradation pathway we consider the two latter cases only, due to the complexity of the involved biochemical reactions. If we assume that the probability of occurrence of each biochemical event depends only on the current state (Markov property), then we can model the stochastic processes describing the life of mRNA molecules as Markov chains.
Thus, following [2,3,4,5], we can write the lifetime probability density ϕ as where S 0 is a square n×n matrix. The off-diagonal elements of S 0 are the transition rates between the transient states, and the diagonal elements are composed of the negative sum over all outgoing transition rates from each transient state. The column vector S A contains all the transition rates governing the transitions from the transient into the absorbing states. The row vector of length n ⃗ e 0 has a one in the position corresponding to the initial condition and zeros elsewhere.
As discussed in Ref. [5], we must take into consideration that at the beginning of the experiment, the expression of mRNA is at the steady state. Thus, not all mRNA molecules are found at the initial state 0; some of them will be in other states because they are older. If the transcription of mRNAs was ongoing long before the beginning of the decay experiment, one can assume that the age distribution of the population of mRNA is at steady state 1 . In [5] it was shown that when the age distribution is at steady state, the fraction of mRNAs remaining ∆t time units after the stop of transcription is given by the function where E(T ) is the average lifetime, namely the average value of T given by and Φ is the probability function of T given by where ϕ(t) can be computed from Eq. (1). An important technical point for the computation of Λ is that in Eq. (2), E(T ) is just the normalization factor in order to have Λ(0) = 1. The networks used in the main article are sufficiently simple, so that the calculation of ϕ can be done easily (see below), without resorting to the relatively complex formulation given in Eq. (1). Nevertheless, Eq. (1) is the general form of the function, which depends on the values of the rates in the matrix S 0 . When Eq. (2) has been written in terms of the underlying rates, the rates can be determined by means of a nonlinear fitting of the log of Λ(∆t) with the log of the data. The use of the logged data is justified by the common assumption that the experimental noise is multiplicative. Eqs. (1) and (2) include the case when there is just one transient state, as a limiting case. This occurs when the decay pattern follows an exponential function whose parameter is the rate of degradation, and also the inverse of the average lifetime.

Fitting functions: General aspects
All our fits were performed for two or three parameters and compared with the exponential fit. Based on the Akaike Information Criterion (AIC) corrected for small sample sizes, we found that the fits with the exponential function did not perform sufficiently well to be considered in this study. Since all our networks are comprised of irreversible transitions only, the computation of Λ can be split in the sub-tasks of computing the contribution of irreversible chains such as the one in figure 1. For each of these chains, the probability density ϕ n for the absorption time in state n can be computed as the convolution of n independent exponential functions, and it is given as if all the ω j are different from one another. In the networks considered in this work we just need the two cases corresponding to n = 1, i.e. the exponential function, and n = 2 because all other cases can be cast in form of a linear combination of these two cases. The cumulative functions associated to the ϕ n are thus given by again under the assumption that all rates are different. As we have seen after fitting the data, this assumption is never violated. Assume now that a given network has, say, three possible paths from the initial state to the absorbing state. Assume also that these paths are made of only irreversible transitions. Let us call these three paths a, b, and c and let us assume that these three paths are characterized each by its own set of ω's and that the probability of each path to be taken is p a , p b and p c , respectively. Let ϕ a , ϕ b and ϕ c be the probability densities conditioned on each of the paths separately. Then, the total probability density ϕ is given by and the cumulative probability function is given by so the decay function Λ from Eq. (2) reads where x = a, b, c, and the proportionality constant must be fixed so that Λ(0) = 1. If we now write the Φ x as Φ x + 1 − 1 and we substitute, we obtain where we have defined The formulation in Eq. (10) turns out to be a quite useful, since our Φ n expressed in Eq. (6) take an explicit relatively simple form once put in Eq. (11). The vector ⃗ ω x is the set of all ω's along the path x. Upon simple explicit integration we have thus, for n = 1 and n = 2 the two functions and to be used in the explicit calculations that follow.

Fitting functions: Specific forms
This section delivers the explicit form of the functions used in the fitting of the data.
• The generic two-state network in Figure 3 in the main article has  Table 1: Experimental data as extracted from the plot published in [1]. The various experimental conditions are given in the first row: (-) = negative control, when the miRNA is not expressed; R = when only miRNA is expressed but both NOT1 and PAN3 are not expressed; RP = data when NOT1 is knocked down; RN = data when PAN3 is knocked down; RNP = all factors are present, no knock downs.
where in both cases ω 1 = λ + µ takes into account that the dwell time on state A is exponential with parameter λ + µ. This function was used to fit the data of the negative control (green line) and therefore to fix the rates λ, µ and ν. These three rates are kept fixed in all other networks.
• The networks in Figures 6, 8, 10 in the main article share the same structure and lead to where λ ′ = λ R , λ RN , λ RN P and µ ′ = µ R , µ RN , µ RN P in Figures 6, 8, 10, respectively. To fit the data, as was done in Figures 7, 9, and 11 in the main article, we have kept the three rates λ, µ and ν fixed as they have been determined in the fit of the negative control. This means that the Λ(∆t) given in Eq. (15) was used to fit the two additional rates λ ′ and µ ′ for each of the following three data sets.

Data used in the study
The experimental data used in this study have been extracted from [1]. We report the numerical values in table 1: the first column is the measurement time expressed in minutes; the second column is the data for the negative control, when the miRNA is not expressed, which was fitted with Eq. (14); the third column is the data when both NOT1 and PAN3 are knocked down (only miRNA is expressed) and was fitted with Eq.  (15); the fourth column is the data when NOT1 is knocked down: since this data is very close to the data in the third column, an additional fit of it was not performed as these two sets of data cannot be distinguished from each other; the fifth column is the data when PAN3 is knocked down and was fitted with Eq. (15); the last column is the positive control, when all factors are expressed and was fitted with Eq. (15).

Fitting procedure and algorithm
Our fitting procedure consisted in finding the set of parameters ⃗ ω that minimized the square deviation of the logarithm of the data, as where "data" is any of the columns in the table above and Λ is either Eq. (14) or Eq. (15) as explained above. The choice of the logarithm is dictated by the assumption that the noise is multiplicative. The minimization was done by using standard MATLAB routines ("lsqnonlin" and "nlinfit") leading to the same results. By using "nlinfit" we were able to estimate also the 95% confidence interval, reported together with the estimated values in the table 2: Since the data concerns averages whereas the biochemical model for the degradation pathway is based on a single molecule perspective, it is to be expected that the confidence intervals are relatively large. In addition, the data used had to be extracted graphically from the plots published in [1] and thus are expected to contain random errors related to reading the data from the plots. Despite these shortcomings, any other of the alternative  Table 3: Steady state properties. For each condition we compute the average lifetime of the mRNA (second column), the percent of molecules that are degraded via miRISC (third column) and the distribution of mRNA among the three main biochemical states of figure 2.
models we have tested were not able to fit the data (not shown), thus indicating that the networks proposed in this study are structurally correct.

Steady state properties and statistics
We can associate to each path of the type drawn in figure 1 its average absorption time given by which allows us to compute the average lifetime associated to each of the networks discussed in our paper. Therefore, the average lifetime associated to the decay pattern described by Eq. (14) is given by where the upper index (-) indicates that this average lifetime refers only to the negative control decay data. Conversely, the average lifetime associated to the decay patterns described by Eq. (15) is given by where the upper index (RX) indicates that this refers to all analyzed networks where miRNA is involved together with none, one or both factors PAN3 and NOT1. When the rates are known one can compute the percent of target mRNA that are involved in miRNA mediated degradation, shown in the third column of  Another quantity of biological interest is the percent of mRNA in the different biochemical states of our network at steady state. The three states of the network are drawn in figure 2 and the stationary probabilities π 0 , π 1 and π 2 can be computed from the Master equation by redirecting the arrows with rates µ, ν and µ ′ towards the initial state 0. This leads to the stationary distribution which have a clear limit when λ ′ → 0 for the negative control case. The values of the probability distribution can be used to estimate how many mRNA should be found in the cell in the three different biochemical states, 0, 1 and 2 in figure 2. Notice that the % of mRNA degraded through miRNA targeting can also be computed in terms of normalized fluxes: consistent with the values computed earlier, given in the third column of table 3.
2 Why the originally hypothesized pathway is not consistent with the experimental data There are several reasons why the biochemical network proposed in [1] is not an adequate model for miRNA-mediated mRNA degradation: 1. The network proposed in [1] does not contain all the state transitions that are  Figure 3: Fit of the data of the decay pattern when PAN3 has been knocked down. The model used here is the one of [1] complemented with a pathways for alternative routes to degradation. In our interpretation this means that λ ′ = λ R in figure 2 is fixed with the data of the experiment with both PAN3 and NOT1 knocked down, while also λ, µ and ν are kept fixed. The only parameter left free to be fixed with this set of data is therefore µ ′ . However, since too few mRNAs are targeted by miRISC alone, the extra parameter µ ′ is not sufficient to fit the data. necessary for fitting both the negative control data and the data where the miRNA is knocked-down.  [1] complemented with a pathway for alternative routes to degradation. In our interpretation this means that λ ′ = λ R in figure 2 is fixed with the data of the experiment with both PAN3 and NOT1 knocked down, while also λ, µ and ν are kept fixed. The only parameter left free to be fixed with this new set of data is therefore µ ′ . However, since too few mRNAs are targeted by miRISC alone, the extra parameter µ ′ turns out to be not sufficient to fit the data.