MIRELLA: a mathematical model explains the effect of microRNA-mediated synthetic genes regulation on intracellular resource allocation

Abstract Competition for intracellular resources, also known as gene expression burden, induces coupling between independently co-expressed genes, a detrimental effect on predictability and reliability of gene circuits in mammalian cells. We recently showed that microRNA (miRNA)-mediated target downregulation correlates with the upregulation of a co-expressed gene, and by exploiting miRNAs-based incoherent-feed-forward loops (iFFLs) we stabilise a gene of interest against burden. Considering these findings, we speculate that miRNA-mediated gene downregulation causes cellular resource redistribution. Despite the extensive use of miRNA in synthetic circuits regulation, this indirect effect was never reported before. Here we developed a synthetic genetic system that embeds miRNA regulation, and a mathematical model, MIRELLA, to unravel the miRNA (MI) RolE on intracellular resource aLLocAtion. We report that the link between miRNA-gene downregulation and independent genes upregulation is a result of the concerted action of ribosome redistribution and ‘queueing-effect’ on the RNA degradation pathway. Taken together, our results provide for the first time insights into the hidden regulatory interaction of miRNA-based synthetic networks, potentially relevant also in endogenous gene regulation. Our observations allow to define rules for complexity- and context-aware design of genetic circuits, in which transgenes co-expression can be modulated by tuning resource availability via number and location of miRNA target sites.


INTRODUCTION
Engineering mammalian cells with synthetic regulatory networks to obtain novel functionalities with predictable behaviour, requires a deep understanding of the dynamic interactions between the genetic circuits and the intracellular context in which they are intended to operate.
We recently showed that when exogenous DNAs are transiently delivered to mammalian cells, they compete for limited shared transcriptional and translational resources, reshaping RNA and protein levels, and leading to the coupling of otherwise independent genes (1,2). This becomes a pervasive problem, either when implementing regulatory circuits, or in attempts to carry out studies using system perturbations (i.e. overexpression or downregulation of a transient gene). The observation that synthetic circuits impose a gene expression burden to their host cells (1,3,4), prompted the development of 'context-aware' gene networks in which incoherent feedforward loops (iFFLs) that use biomolecular controllers such as endoribonucleases (2), or endogenous and synthetic miRNAs (1) were successfully exploited as burden mitigators. In miRNA-based iFFL, we observed that the protein levels of two independently expressed genes, one of which was regulated by miRNA, were strongly linked to the miRNA activity (1). As a point in case, using two co-expressed fluorescent proteins, EGFP and mKate, with mKate levels linked by design to the endogenous miRNA-31 (miR-31), we showed that the higher the number of miRNA target sites (TS) in the mKate UTRs, the stronger its downregulation, and the higher the levels of EGFP (Figure 1A). Conversely, by inhibiting miR-31, mKate levels increased, counterbalanced by reduced EGFP levels (1), supporting the miRNA-dependency of the observed effect. Results were robust to changing plasmid design (cotransfection vs single plasmid) and cellular context (1).
Given the biological importance and high applicability of miRNAs to synthetic circuits, we sought to investigate the molecular mechanisms involved in miRNA-dependent resource distribution. We postulate that this understanding will enable more precise circuits' design with enhanced robustness and predictability, and may shed light on the secondary regulatory effect in endogenous pathways.
miRNAs are small non-coding RNAs produced from transcripts with stem-loop structures which undergo processing both in the nucleus and cytoplasm to be converted into mature, 21-26 nucleotides-long miRNAs. Mature miR-NAs are assembled into the RNA-induced silencing complex (RISC) and bind their target mRNAs by base pairing usually to their 3'UTR or 5'UTR (5,6). Upon binding, miRNAs modulate their target through mRNA degradation and/or translational repression (7,8). Target sites (TS) can be fully or partially complementary to the miRNA. In the former case, the target is degraded through endonucleolytic cleavage (9), while in the latter translational repression dominates and transcript degradation occurs after deadenylation (7,8). Typically, in mammalian cells miRNAs bind to the 3'UTR of the endogenous target mRNA, and have non-perfect complementarity to the TS (7). Thus, their main mechanism of action is to repress translation. Moreover, miRNAs are often found in endogenous feedforward or negative feedback loops exploiting additional functions such as buffering gene expression against noise or fluctuations in external inducer concentration (10)(11)(12)(13).
In synthetic biology, miRNAs have been repurposed as a versatile tool to build cell-specific devices, and have been largely used to create cell classifiers with biotechnological or biomedical applications (14,15) or to modulate the expression of the genetic devices (16,17). In these applications, to achieve strong downmodulation of the target genes, and increase the sensitivity to small concentrations of miR-NAs, perfectly complementary target sites are typically used (18).
Here, we use a two-gene reporter system (EGFP hereafter named capacity monitor, and miR31TS-mKate, hereafter named miTarget), along with a mathematical model (MIRELLA) that qualitatively captures posttranscriptional events, to explain the effect of miRNA (MI) REgulation on intracellular resource aLLocAtion, effectively identifying key processes responsible for miRNAbased burden mitigation in mammalian cells ( Figure 1A, B).
MIRELLA builds on an existing modelling framework (1) and considers that mRNA translation and degradation use pools of shared resources, among which we account as main players ribosomes and RNases, respectively. MIRELLA replaces reaction rates that involve shared cellular resources with effective reaction rates that account for the availability of each individual resource pool according to the overall gene expression demand ( Figure 1B). We use MIRELLA to predict the effect of miRNA regulation on resource availability considering the strength of the downregulation (number of target sites), and the effect of TS location (i.e. as part of the 5'UTR or 3'UTR). The model suggests that miRNA-mediated downregulation of the target mR-NAs causes a redistribution of translational resources (i.e. ribosomes) and impacts the RNA degradation machinery, overall contributing to a change in protein expression levels ( Figure 1C). We then experimentally validated the model predictions for the expression of exogenous and endogenous genes.
Supported by the synergistic use of a mathematical model and experiments, our findings contribute to a deeper understanding of the mechanisms of miRNA operations in synthetic networks, which enables the resource-aware design of genetic circuits. Moreover, our results provide insights into secondary effects of miRNA regulation that might be potentially relevant also in endogenous gene regulation.

DNA cloning and plasmid construction
Plasmid vectors carrying gene cassettes were created using In-Fusion HD cloning kit (Clontech), or digestion and ligation. Reaction included 1:2 molar ratio of plasmid backbone:gene inserts starting with 100 ng of vector backbone digested with selected restriction enzymes. All plasmids used in this study consist of a constitutive promoter driving the gene of interest. All plasmids used in this study were confirmed by sequencing analysis and are listed in Supplementary Table 2.

Cell culture
H1299 cells were maintained in Roswell Park Memorial Institute medium (RPMI, Gibco) supplemented with 10% FBS (Atlanta BIO), 1% penicillin/streptomycin/Lglutamine (Sigma-Aldrich) and 1% non-essential amino acids (HyClone). The cells were maintained at 37 • C and 5% CO 2 . ) that capture the availability of that resources according to the overall demand from competing genes (modelled via the general function f (. . .)). (C) Effect of miRNA activity on protein expression in a finite-resource context. Control (noTS): protein expression in the absence of miRNA regulation. TS in the 5' or 3'UTR: the slicer activity of miRNA-RISC complex triggers mRNA degradation of the miTarget (red) causing a queueing effect on the degradation of other mRNAs that results in capacity monitor (green) protein accumulation. TS in the 5'UTR: in addition to slicer activity, miRNA binding inhibits translation initiation, freeing up translational resources to the benefit of other transcripts (e.g. the capacity monitor). The resulting effect is stronger downregulation of miTarget and higher capacity monitor levels when TS are placed in the 5'UTR of the miTarget.

Transfection
Transfections were carried out in 24-well plate format for flow cytometry analysis, in 12-well plate format for flow cytometry and qPCR analyses run on the same biological replicate or in 10 cm dishes for polysome profiling. H1299 cells were transfected with Lipofectamine® 3000 (Thermo Fisher Scientific) according to the manufacturer's instructions and 300 ng of total DNA for each sample for a transfection in 24-well plates and scaled up for larger formats. Details on transfections are provided in Supplementary  Table 1.

Flow cytometry and data analysis
Cells were analysed using a BD FACSAria™ cell analyser (BD Biosciences) using 488 nm and 561 nm lasers. Cells transfected in 12-well plates were washed with DPBS, detached with 100 l of trypsin-EDTA (0.25%) phenol red and resuspended in 600 l of DPBS (Thermo Fisher Scientific). 200 l of cell suspension were used for flow cytometry and 400 l for RNA extraction. For each analysis, 10 000+ events from each sample were recorded and data were normalised with three compensation controls: unstained (wild-type cells), and single colour controls (mKate only, EGFP only). Fluorescence intensity in arbitrary units (AU) was used as a measure of protein expression. Population of live cells was selected according to FCS/SSC parameters. Data analysis was performed with Cytoflow. For each sample, we gated the population of live cells and then the EGFP + mKate + cells (Q2 quadrant in Supplementary Figure 9b). Within this population we calculated the geometric mean (Geo-Mean) of mKate and EGFP.

Polysome profiling
Polysome profiling was performed following the protocols described in (19,20). To obtain the cytoplasmic lysates, cells were treated with cycloheximide (10 g ml −1 ) for 3-4 min and then lysed in 300 l of cold hypotonic lysis buffer (19). To remove nuclei, mitochondria and cellular debris, the lysates were centrifuged at 4 • C for 5 min at 20 000 g. To separate ribosomal subunits, ribosomes and polysomes from other cytoplasmic molecules, the supernatant was loaded on a 10-40% (w/v) sucrose gradient and centrifuged for 1 h 30 min at 260 000 g at 4 • C in a SW41 rotor using a Beckman Optima LE-80 Ultracentrifuge. Twelve 1 ml fractions were collected and the absorbance at 254 nm was monitored with the UA-6 UV/VIS detector (Teledyne Isco). RNA was purified fraction by fraction using the phenol/chloroform extraction method described in (21). The retro-transcription reaction was performed using the same volume of RNA for all polysomal fractions. The cosedimentation profile of mRNAs was obtained by calculating the percentage (or fraction) of mRNAs in each fraction by qPCR as described in (22).

mRNA extraction and reverse transcription
RNA extraction was performed with E.Z.N.A.® Total RNA Kit I (Omega Bio-tek). Protocol was followed according to the manufacturer's instructions and RNA was eluted in 30 l of RNase-free water to maximise the yield. RNA samples were conserved at -80 • C. The protocol was performed exclusively with RNase free water in an RNase-free environment.
PrimeScript RT Reagent Kit with gDNA Eraser--Perfect Real Time (Takara) was used according to manufacturer's instructions. The protocol was performed on ice in a RNase-free environment to avoid RNA degradation. A negative control without PrimeScript RT Enzyme Mix I was always prepared to be sure that samples were not contaminated with genomic DNA.

qPCR
Fast SYBR Green Master Mix (Thermo Fisher Scientific) was used to perform qPCR of cDNAs obtained from 500 ng of RNA and diluted 1:5. Samples were loaded in Mi-croAmp™ Fast Optical 96-Well Reaction Plate (0.1 ml) and the experiment was carried out with a 7900HT™ Fast machine. Each well contained 10 l of final volume (5 l SYBR Green Master Mix 2X, 2 l ddH2O, 1 l of each primer, 1 l of template). Also, a control without template (blank) was set. Primers were designed to amplify a region of 60-200 bp (Supplementary Table 3) and with a temperature of annealing between 50 • C and 65 • C. Data were analysed using the Comparative Ct Method according to Applied Biosystems Protocols.

Modelling
We constructed a deterministic ODE model to qualitatively capture post-transcriptional events in order to identify key processes responsible for miRNA-based resource reallocation. Our ODE model is based on a previously published resource-aware modelling framework (1). Full details of its formulation and parameterisation can be found in Supple- mentary

Deterministic simulations of the ODE model
All the simulations of the ODE model were run using Python 3 (v. 3.8.13). The scipy.integrate.solve ivp function from the SciPy library (v. 1.8.0) was used to numerically integrate the ODE model. More specifically, we used the Radau method (stiff ODE solver) with an absolute tolerance of 10 −9 , and a relative tolerance of 10 −6 . To simulate the ODE model at steady state, the time span was set to [0, 10 000] h. Hence, steady state was taken as the value for each molecular species at the end of a numerical simulation of 10 000 h.
The parameter values used to simulate the ODE model were chosen as described in Supplementary Note 6 and reported in Supplementary Table 5 (for simulations of the  translational resource reallocation) and Supplementary Table 7 (for simulations of the degradation resource reallocation). All plots were generated in Python using the Matplotlib library (v. 3.5.1). The code to run all the simulations and the model is available publicly as Jupyter notebooks (https://github.com/giansimone/MIRELLA/).

Analytical characterisation of the translational resource reallocation
To quantitatively characterise the reallocation of the translational resources at steady state, we derived an analytical solution of the ODE model as reported in Supplementary Note 1. Full details of the analytical solution derivation can be found in Supplementary Note 4. The steady-state expression levels for both the miTarget (p T ) and the capacity monitor (p C ) are as follows: where ρ T is the resource demand coefficient for the miTarget, ρ Q T is the resource demand coefficient for the miTarget:miRNA complex, and ρ C is the resource demand coefficient for the capacity monitor. The boolean parameter σ ∈ {0, 1} captures the location of the miR-TS at either the 5' (σ = 0) or 3' (σ = 1) UTRs. The analytical expressions of the resource demand coefficients are as follows: where κ T and κ C are the effective dissociation constant for the miTarget and capacity monitor, respectively (see Supplementary Note 1 for their definitions). All other parameters are described in Supplementary Table 5.

Model fitting
The model fitting shown in Supplementary Figure  3 was performed using Python 3 (v. 3.8.13). The scipy.optimize.differential evolution function from the SciPy library (v. 1.8.0) was used to find the model parameters that best fit the predicted steady-state expression levels for both the miTarget and capacity monitor to the experimental data reported in Figure 2B. The error dis-tance between the predicted and the experimental data was evaluated using the following loss function: is the vector that contains the different mean values for the miTarget (mKate) and the capacity monitor (EGFP) in each condition as reported in Figure 2b, T is the vector that contains the predicted values in the different simulated conditions, and ϑ is the vector that contains the model parameters that have to be identified via the model fitting.
An L 2 -regularisation term was added to the loss function to prevent an ill-conditioned parameter estimation problem.
To keep all the model parameters within the same order of magnitude, the regularisation hyperparameter was set to λ = 0.001. The plot in Supplementary Fig. 3 was generated in Python using the Matplotlib library (v. 3.5.1). The code to perform the model fitting is available publicly as a Jupyter notebook (https://github.com/giansimone/MIRELLA/).

A modelling framework explains the relation between gene downregulation by microRNAs and cellular resources reallocation
We previously observed that miRNA-based iFFL circuit designs can be used to reduce the resource-based coupling of two co-expressed genes, and that miRNA-driven gene downregulation is associated with increased expression of co-encoded, independent genes (1). To understand the mechanisms underlying this indirect effect on independent genes expression, we developed MIRELLA, a deterministic resource-aware model that captures the resourceconstrained co-expression of two constitutive genes when one of them is downregulated by an endogenous miRNA (Figure 2A). With a focus on miRNA activity as a key player in iFFL-based burden mitigation, MIRELLA extends existing models of biochemical reactions by taking into account the effects of shared cellular resources. This is achieved through the use of effective reaction rates that depend on cellular resources demand, like previously done in (1), while explicitly accounting for ribosomes and RNases resource pools ( Figure 1B). In what follows, we provide a brief overview of the model, whilst its full details can be found in Supplementary Note 1.
Since miRNAs act post-transcriptionally, the model focuses on the processes that contribute to proteins' expression, namely miRNA-dependent mRNA downregulation, mRNA degradation and mRNA translation. The driving hypothesis is that the observed increased expression of an independent gene following the downregulation of the miRNA target is likely due to the reallocation of shared cellular resources involved in post-transcriptional processes. The model thus consists of a set of ordinary differential equations (ODEs) that capture the relative gene expression levels of a constitutive gene which is modulated by an endogenous miRNA (miTarget-- Figure 2A), and a second constitutive gene, whose expression levels reflect variations in the availability of shared cellular resources (capacity Nucleic Acids Research, 2023, Vol. 51, No. 7 3457  Table 4, whilst all the model parameters--including the numerical values used for the presented simulation results--are summarised in Supplementary Table 5. monitor-- Figure 2A). The miTarget is a mRNA encoding for a red fluorescent protein (mKate) that includes target sites for miR-31 in its 5' or 3' UTR, whereas the capacity monitor is a constitutively expressed green fluorescent reporter (EGFP). For each of these two genes, the model captures the essential features of transcription, translation, degradation, and also the interactions between genes and shared cellular resource pools as illustrated in Figure 2A. The model assumes that the shared cellular resource pool for translation (here comprising only ribosomes for simplicity) is constant, whereas the pool for RNA degradation (e.g. RNases) is initially considered unlimited. This enables us to initially neglect the demand for RNA degradation resources and to focus on how miRNA-driven regulation impacts ribosomes reallocation (Figure 2A and 3A,  B). In a second step, we consider a finite pool of RNases, which are reallocated depending on degradation demand from the co-expressed genes ( Figure 4A 1). These rates change dynamically according to the translational resource demand, which depends on the mRNA expression levels m T and m C , and on the effective dissociation constants κ T and κ C (see Supplementary Note 1 for further details). The model considers two main modes of action of miRNAs on their target genes, namely regulation of translation initiation (23,24) and mRNA degradation (25). Since synthetic circuits commonly use perfectly complementary TS to maximise the downregulation, we consider mRNA degradation the main outcome of miRNAs activity. However, when located in proximity of the AUG at the 5' UTR, we additionally consider a steric hindrance effect of miRISC complex that competes with ribosome binding and impairs mRNA translation (26) ( Figure 1C).

Mass-action kinetics
To capture the binding interactions between the miTarget mRNA and its cognate miRNA, the model considers a further molecular binding reaction as illustrated in Figure 2A. The strength of the miRNA regulation is modelled via the characteristic parameter associated with the binding reaction between the miTarget mRNA and its cognate miRNA, that is the miRNA binding constant η + ( Figure  2A). η + can be used as a proxy for capturing the effect of different repetitions of TS in the 5' or 3' UTRs of the target mRNA. Hereafter, the binding constant η + is assumed to be proportional to the number of TS (i.e. we do not consider cooperation effects in miRNA-driven regulation). The '3' UTR model' (shaded orange region in Figure 2aA) considers miRNA regulation via the mRNA degradation rate β Q T , whilst the '5' UTR model' (shaded blue region in Figure 2A) considers miRNA regulation via both the mRNA degradation rate β Q T and the effective translation rate γ E f f T, Q (Supplementary Note 1). We assume first-order miRNA-mediated degradation of the target mRNA miTarget, which results in an increase of the mRNA degradation rate by a factor λ Q T > 1, that is β Q T = λ Q T × β T . When the TS are within the 5' UTR, the model assumes that the effective translation rate γ E f f T, Q is equal to zero since, in this case, miTarget:miRNA transcripts cannot be translated due to steric hindrance.
To check if our model correctly predicts protein expression upon miRNA regulation, we ran simulations of the resource-aware model and compared the results with the experimental data published in (1). Data represent the protein levels of miTarget and capacity monitor in H1299 cells, which naturally express high levels of miR-31. miTarget and capacity monitor are under the regulation of a bidirectional constitutive CMV promoter (Supplementary Figure 9a). The miTarget includes either 1 or 3 TS in the 3' or in the 5' UTR. The higher the number of TS, the stronger the repression exerted by miR-31, and consequently, the higher the capacity monitor levels ( Figure 2B).
We simulated the presence of multiple miRNA TS in both the 3' UTR and 5' UTR by considering a range of reasonable characteristic values for the miRNA binding constant η + (Supplementary Table 5 and Supplementary Note 6), and monitored at the steady-state the values of the molecular species of the system (Materials and Methods). Simulation results confirmed that the resource-aware model can recapitulate the steady-state protein levels of the miTarget and the capacity monitor observed in the experimental data ( Figure 2C). We found that miRNA binding constant η + is positively correlated with the steadystate protein levels of the capacity monitor. As expected, the simulation results showed that when miRNA TS are located in the 5' UTR there is stronger downregulation of the miTarget followed by higher capacity monitor levels ( Figure 2C).
We next solved the model analytically to obtain insights into resource demands that account for the steady-state levels of both the miTarget and the capacity monitor proteins (Eq. (1), Materials and Methods; full details in Supplementary Note 4). We found that these levels depend on the resource demand coefficients ρ T (for miTarget translation, Eq. (2a)), ρ C (for capacity monitor translation, Eq. (2c)), and ρ Q T (for miTarget:miRNA translation, Eq. (2b)). We noticed that the resource demand coefficient ρ Q T is nonzero only when TS are present in the 3'UTR (Eq. (2b)). We also observed that the miRNA regulation directly alters the resource demand coefficients associated with mi-Target translation, although in different ways. More specifically, the coefficient ρ T (Eq. (2a)) is downregulated by both the miRNA binding constant η + and the miRNAenhanced degradation rate β Q T . In contrast, the demand coefficient ρ Q T (Eq. (2b)) is upregulated by the miRNA binding constant η + and downregulated by the miRNA-enhanced degradation rate β Q T . Therefore, the miRNA-driven miTarget degradation reduces the resource demand coefficients, and hence the amount of miTarget mRNA that is translated into miTarget proteins. This effect leads to a redistribution of the translational resources on the capacity monitor transcripts, which in turn increases the amount

Regulation by miRNA binding to the 5'UTR results in ribosomes reallocation and modified translational profiles
To explore the effect of miRNA-driven reallocation of translational resources we first used our resource-aware model to predict the amount of free and translating ribosomes at steady state. Based on the results of the simulation reported in Figure 2C, the pool of free ribosomes positively correlates with the miRNA binding constant η + Nucleic Acids Research, 2023, Vol. 51, No. 7 3461 ( Figure 3A), and with the location of miRNA-TS. Furthermore, the ribosomal density defined here as the number of translating ribosomes per total number of transcripts, positively correlates with the miRNA binding constant η + for both the miTarget and the capacity monitor ( Figure 3B). This suggests that the miRNA-driven downregulation of the miTarget frees up translational resources, which can be deployed to the translation of other transcripts, including co-expressed ones such as the capacity monitor. We investigated what may cause the ribosomes reallocation by analytically solving the model at the steady state (Supplementary Note 4). We derived the steady-state level for the pool of free ribosomes (r ) and found that it is inversely proportional to the sum of the resource demand coefficients, that isr ∼ Note 4). Since miRNA regulation does not affect the capacity monitor resource demand coefficient, a reduction in the total miTarget resource demand (i.e. ρ T + ρ Q T ) corresponds to an increase in the amount of free ribosomes (Supplementary Note 4). We then analytically calculated the steady-state ribosomal densities for the capacity monitor transcripts and found that these correlate with the amount of free ribosomes (Supplementary Note 4). Of note, our simulations and the analytical solution show an increased number of ribosomes per transcript also for the miTarget mRNA when the miR-TS are located at the 5' UTR ( Figure 3B, Supplementary Note 4). Even if counterintuitive, this suggests that the increased amount of available ribosomes is to the benefit of miTarget transcripts that did not undergo miRNA-mediated downregulation.
Next, we experimentally validated the model predictions. We performed polysome profiling in H1299 cells transfected with miTarget and capacity monitor to observe the changes in co-sedimentation of the reporter mRNAs with polysomes (19,20) following miRNA modulation. We considered several conditions, including no miRNA regulation and TS inserted in the 5' or 3' UTR ( Figure 3C, D). After sucrose gradient separation of cytoplasmic lysates, it is possible to follow the sedimentation of transcripts with free cytosolic light components (ribonucleoproteins, RNPs), ribosomal subunits (40S and 60S) and monosomes (80S) -all associated with non-translating particles -and with polysomes comprising translating transcripts bound by multiple ribosomes (Figure 3C, D). In line with model predictions, both the capacity monitor (EGFP, Figure 3C, bottom) and the miTarget (mKate, Figure 3C, top) exhibit modified co-sedimentation profiles upon miR-31 modulation as compared to the control (no miRNA regulation, grey line). In addition, we monitored the co-sedimentation profiles of endogenous transcripts, i.e. eIF4E, CCNA2 and GAPDH, that we previously observed to be impacted by resource competition. Following miRNA regulation of mi-Target, we found that they also benefit from ribosome reallocation ( Figure 3D). Our data indicate that for all analysed genes, there is a shift towards the polysomal fractions (peaks 9-10) that is more pronounced when the miR-31 TS are placed at the 5' UTR ( Figure 3C, D, left) as compared to the 3' UTR ( Figure 3C, D, right) of the miTarget. These results are consistent with our hypothesis stating that, due to the physical proximity of the TS to the AUG, the miRISC complex interferes with the binding of ribosomes to the target mRNA, resulting in their re-allocation on other transcripts.

microRNA activity induces a queuing effect on mRNA degradation
Along with ribosome redistribution, we used MIRELLA to investigate the impact of miRNA-driven regulation on the RNA degradation machinery. To this end, we modified our model to include a finite pool of RNases (Supplementary Note 5), recapitulating the essential features of transcription, degradation, and interactions between transcripts and the degrading resource pool (RNases) as illustrated in Figure 4A. We then simulated the degradation of an endogenous transcript upon halting of transcription (i.e. all transcription rate constants were set to zero after the system reached its steady state) for different reasonable values of the miRNA binding constant η + (Supplementary Table 7 and Supplementary Note 6).
The competition for the shared pool of RNases is captured via the effective mRNA degradation rates derived for the miTarget (β  Supplementary Fig. 2 with simpler first-order decay reactions ( Figure 4A and Supplementary Note 5). These rates change dynamically and depend both on the expression levels of the degrading mRNA:RNase complexes (s T and s E , Supplementary Fig. 2) and the association and dissociation constants between the RNase and the mRNA strands (Supplementary Note 5 for further details). We assume that the miRNA-driven decay increases the rate of binding between the miTarget transcript and the molecular species of the degradation machinery by increasing the association constant β + T to the new value β + T, Q (Supplementary Fig. 2). This effect produces an increase in the effective mRNA degradation rate β  Figure 4B, extended data presented in Supplementary Fig. 4) and showed that the endogenous degradation dynamics are slower when compared to the control (η + = 0). Moreover, the endogenous degradation profile presents two different decay phases when the miRNA binding constant η + is greater than ∼2 nM −1 h −1 ( Figure 4B and Supplementary Figure 4b). Specifically, the endogenous mRNA levels remain almost stable during the first decay phase, whereas they degrade quicker during the second phase. We hypothesise that the increased degradation of the target by miR-31 ( Supplementary Figure 4a) sequesters a significant portion of the degradation machinery ( Supplementary Fig. 4c, d), which is thus not available for the degradation of non-target mRNAs (Figure 4A and Supplementary Figure 4b, e). We also observed that the decay dynamics of the non-targeted transcripts become faster once the majority of miTarget transcripts are degraded ( Supplementary Figures 4 and 5), probably due to the reallocation of degrading resources from the mi-Target mRNAs to the non-target mRNAs ( Supplementary  Figure 4d, e).
To experimentally validate the model predictions, we treated cells with 5,6-dichloro-1-beta-D-ribobenzimidazole (DRB), a transcription inhibitor (27) and measured the mRNA half-life of the endogenous genes CCNA2 and eIF4E, whose half-life are reported within the range 2-4 h (28) (Figure 4C, D). We selected these genes since their mRNA half-life is much shorter than that of the capacity monitor EGFP, which is roughly >8 h (29) and thus does not permit to capture faster decay dynamics.
We transfected H1299 cells with mKate encoding plasmid with noTS or 3 miR-31 TS at either the 3' or the 5'UTR along with a constitutively expressed EGFP. Over a timeframe of 4 h, CCNA2 and eIF4E mRNAs exhibit a longer decay time when mKate is flanked by miR-31 TS as compared to the control (no miR-31 TS) ( Figure 4C, D). More specifically, the two mRNA species decrease immediately after treatment only in the control, while in the case of 3TS at the 3' and 5' UTRs, they remain almost stable for 1 h 30 min after treatment, undergoing a significant decrease after 2-3 h. The qualitative trend of the 3TS 3' and 3TS 5' samples are very similar and in agreement with our hypothesis that miRNAs enhance target degradation at the same rate when they bind to 3'-or 5'-UTR TS. Finally, the experimental measurements show that the two mRNA species reach a plateau 4 h post-DRB treatment, in agreement with the model predictions. Although the predicted and experimental decays exhibit the same dynamics, we find that the endogenous mRNAs are not completely degraded upon treatment with DRB, which we argue is a limitation of the mRNA half-life assay. In fact, we measured CCNA2 and eIF4E mRNAs up to 6 h post-DRB treatment in wildtype cells and observed that it reaches about 50% of the initial concentration for both genes (Supplementary Figure 8). Overall, our results confirm the model predictions and are in agreement with our hypothesis of a prominent queuing effect on the mRNA degradation pathway as a consequence of the increased degradation of the miTarget by miR-31.

DISCUSSION
miRNAs are fundamental building blocks of posttranscriptional control of gene expression and help to finely control genetic circuits and precisely regulate endogenous pathways. We observed that when a synthetic transcript is downregulated by an endogenous miRNA, a synthetic co-expressed mRNA is indirectly upregulated. Here, we combined mathematical modelling and experimental analysis to understand important mechanisms underlying the effect of miRNAs on resource allocation in mammalian cells, when their regulation is embedded in synthetic circuits. This understanding is instrumental to generate more detailed guidelines for an informed and rational design of gene circuits, and may also provide insights into endogenous gene regulation.
We had previously observed that exogenous genes delivered in mammalian cells compete for intracellular resources. This effect was observed, albeit to a different extent, regardless of genes amount and cell types (1).
Since miRNAs function at post-transcriptional level, we developed a mathematical model, MIRELLA, that considers effective reaction rates of biological processes to account for the availability of post-transcriptional shared resources. Specifically, the model considers one pool of mRNA translation resources and one pool of mRNA degradation resources, for simplicity represented by ribosomes and RNases, respectively (Figures 2A and 4A). Although the aforementioned biological processes are regulated by many molecular species, the key cause of resource competition in mammalian cells still remains an open question. Here we focused on ribosomes and RNases for their relevance in these processes and--relatively to ribosomes--for the possibility to experimentally measure them. This allows us to keep the model as simple as possible to study the role of two prime contributors--that is ribosomes and RNases--in gene expression burden. We anticipate that other molecular species could contribute to shaping gene expression as a result of resource competition, such as tRNAs for rare codons in mRNA translation. In this respect, our resource-aware modelling framework can be extended to capture more complex scenarios, for example, by including the resource competition for shared tRNAs in the model's equations. Moreover, our context-aware theory paves the way towards a holistic understanding of cell biology (30) instrumental to engineer more sophisticated biomolecular circuits that can maximise the global system's output (e.g. a protein or a compound of interest) by operating synergistically with the cell physiology (31) The experimental setting is based on two fluorescent reporters driven by constitutive promoters, one regulated by an endogenous miRNA (miTarget) and the other lacking regulation, as a proxy for resources availability (capacity monitor). This design is simple and is meant to quickly appreciate differences in protein expression levels. The modelling framework suggested that miRNA regulation has two major consequences on resource allocation. One is a queueing effect on the RNA degradation pathway as the result of the strong miRNA-mediated slicing of miTarget transcripts. The other is an increased amount of ribosomes becoming available for translating other transcripts.
Synthetic networks that embed miRNA regulation typically employ perfectly complementary TS to maximise the fold change expression of the target gene (18), and it has been generally observed that the abundance of miRNAs with respect to their target affects genetic circuits functionality (32). It was previously proposed that a significant mRNA degradation by miRNA may lead to a 'queueing effect' for the degradation of other mRNAs, decreasing the effective mRNA decay rate (33). We validated this hypothesis by quantifying the mRNA decay of two endogenous genes, CCNA2 and eIF4E, which were previously shown to be impacted by the burden imposed by the expression of synthetic circuits (1). We observed that upon transcriptioninhibition treatment with DRB, CCNA2 and eIF4E are degraded at a slower pace in the samples expressing miTarget as compared to the control lacking the TS. These results suggest that an accumulation of other mRNAs contributes to increased cognate protein levels. This finding may be useful for circuit design considerations, particularly when such systems are devised to study endogenous processes or pathways.
Our experimental findings suggest that ribosomes reallocate upon miRNA regulation and that this effect is more