Modelling Negative Feedback Networks for Activating Transcription Factor 3 Predicts a Dominant Role for miRNAs in Immediate Early Gene Regulation

Activating transcription factor 3 (Atf3) is rapidly and transiently upregulated in numerous systems, and is associated with various disease states. Atf3 is required for negative feedback regulation of other genes, but is itself subject to negative feedback regulation possibly by autorepression. In cardiomyocytes, Atf3 and Egr1 mRNAs are upregulated via ERK1/2 signalling and Atf3 suppresses Egr1 expression. We previously developed a mathematical model for the Atf3-Egr1 system. Here, we adjusted and extended the model to explore mechanisms of Atf3 feedback regulation. Introduction of an autorepressive loop for Atf3 tuned down its expression and inhibition of Egr1 was lost, demonstrating that negative feedback regulation of Atf3 by Atf3 itself is implausible in this context. Experimentally, signals downstream from ERK1/2 suppress Atf3 expression. Mathematical modelling indicated that this cannot occur by phosphorylation of pre-existing inhibitory transcriptional regulators because the time delay is too short. De novo synthesis of an inhibitory transcription factor (ITF) with a high affinity for the Atf3 promoter could suppress Atf3 expression, but (as with the Atf3 autorepression loop) inhibition of Egr1 was lost. Developing the model to include newly-synthesised miRNAs very efficiently terminated Atf3 protein expression and, with a 4-fold increase in the rate of degradation of mRNA from the mRNA/miRNA complex, profiles for Atf3 mRNA, Atf3 protein and Egr1 mRNA approximated to the experimental data. Combining the ITF model with that of the miRNA did not improve the profiles suggesting that miRNAs are likely to play a dominant role in switching off Atf3 expression post-induction.

. Arrows indicate protein activation/upregulation, mRNA transcription or protein translation whilst negative regulation is shown by horizontal lines. The symbol "ϕ" indicates degradation of the respective mRNA/protein/complex as indicated.

S1 Revised original model: Atf3 protein inhibits transcription of Egr1
In the original model in which Atf3 regulates Egr1 mRNA (S6) the action of phosphorylation of the transcription factors driving Egr1 and Atf3 mRNA expression levels were implicitly assumed. In this model revision we explicitly assume the binding of transcription factors to the associated Egr1 and Atf3 genes and their subsequent phosphorylation by phosphorylated ERK. This allows us to include a description of phosphorylated ERK de-phosphorylation which was not considered in our earlier work.

S1.1 Reaction equations
The phosphorylation of MKK by ET-1 is denoted by which subsequently phosphorylates the unphosphorylated ERK MKK-P + ERK k 2 → MKK-P + ERK-P.
(2) ERK-P subsequently de-phosphorylates such that ERK-P is now free to phosphorylate transcription factors of Egr1 (TF Egr1 ) and Atf3 DNA (TF Atf3 ), which we assume are already bound to the DNA of each species, such that and TF Atf3 · DNA Atf3 + ERK-P k 5 k −5 TF-P Atf3 · DNA Atf3 + ERK-P, which are subsequently transcribed to produce the associated mRNAs and TF-P Atf3 · DNA Atf3 Both the phosphorylated bound transcription factors dephosphorylate in time and the mRNAs subsequently degrade mRNA Egr1 Here · denotes a complex and ϕ the degraded entity.
The suppression of Egr1 mRNA transcription by Atf3 is described by Finally the translation of Atf3 mRNA to Atf3 protein and subsequent degradation of the protein are denoted by mRNA Atf3 respectively, where ϕ P denotes degraded protein.
In this work we do not explicitly account for the degraded mRNAs, Atf3 protein or dephosphorylated ERK-P.

S1.2 Mathematical Model
Applying the Law of Mass Action (S8) to equations (1)-(12) yields where e T represents the concentration of ET-1 (denoted e T =[ET-1] and assumed con- Here all non-zero initial conditions are given by

Conservation of T E and T A means
and which along with solving for m P as detailed in (S6) gives the reduced system of equations with the initial conditions

S1.2.1 Parameterisation, model solutions and sensitivity analysis
Model parameters are summarised in Table S1 in Text S1. The revision of our original model (S6) requires knowledge of 8 further parameters. These are: (i) the forward and reverse rates of bound transcription phosphorylation for Egr1 DNA and Atf3 DNA (k 3 , k −3 , k 5 , k −5 ), respectively; (ii) the reverse rate of Atf3 binding to the phosphorylated transcription factor on Egr1 (this was excluded in our original model for reasons of simplification); (iii) the rate of ERK-P de-phosphorylation (d 9 ); and (iv) de-phosphorylation rates for the phosphorylated transcription factors bound to Egr1 and Atf3.
Given the qualitatitve nature of the data available to fit these new parameters ( Figures  1D, 1F and 1G) a fit-by-eye was determined as the best method for fitting the mathematical model to the experimental data. Having determined these values a local sensitivity analysis, whereby each parameter was varied 2 orders of magnitude above and below their initially estimated value (see Table S1 in Text S1), was undertaken in order to evaluate how good a prediction each parameter estimate was. This was likewise undertaken for all the remaining models.
The forward and reverse rates of Egr1 and Atf3 bound transcription factor by ERK-P (k 3 and k 5 ) were obtained by fitting T EP and T AP to the data in Figures 1F and 1G,  respectively. Estimates were then made of the reverse rates (k −3 and k −5 ) as shown in Table S1 in Text S1. ERK-P de-phosphorylation (d 9 ) was obtained by making a good fit-by-eye between the data shown in Figure 1D. An initial estimate of 1.00 ×10 −3 s −1 was made for k −7 , the rate of reverse association of Atf3 for Egr1 DNA, however, this value greatly decreased the levels of Egr1 mRNA below those experimentally observed (see Figure 1B). Utilising our sensitivity analysis we determined a value of 1.00 ×10 −4 s −1 gave a much better fit to the model. The dephosphorylation rates of the Egr1 and Atf3 bound transcription factors (TF-P Egr1 and TF-P Atf3 ) were determined by fitting the model with data shown in Figures 1F and 1G. Finally, a best fit for the Egr1 data was obtained by increasing the rate of both Atf3 and Egr1 transcription 5-fold in comparison to that originally stated in (S6). This increase is feasible given the earlier work used the lower bounded estimate of transcription for both genes.
We also chose to investigate variation of all the remaining parameter values, each of which had been informed from experimental data within the literature as detailed in Table S1 in Text S1. These results revealed no unexpected behaviour. For instance decreasing the rate of Egr1 transcription meant not enough Egr1 mRNA was produced to be discernable; increasing Egr1 transcription did not affect the initially observed results. Decreasing the transcription of Atf3 reduces the amount of Atf3 produced, which in turn increases the amount of Egr1 mRNA (via the negative feedback relationship) such that it no longer fits our experimental data. Decreasing and increasing the rate of Atf3 translation had a similar effect. Likewise increasing/decreasing Atf3 mRNA degradation and likewise Atf3 led to more/less Egr1 mRNA being produced.
Solutions to the revised model are shown in Figure S2 in Text S1.

S2 Model Extension 1 -Atf3 protein inhibits Atf3 transcription
Here we consider Pathway 1 in Figure S1 in Text S1 where Atf3 protein binds to Atf3 DNA to inhibit its own transcription which is in competition with the binding of the transcription factor TF Atf3 which is required for Atf3 transcription TF Atf3 + DNA Atf3 and which is subsequently phosphorylated by ERK-P TF Atf3 · DNA Atf3 + ERK-P Parameter Definition Atf3 mRNA transcription rate. Rate of ERK-P de-phosphorylation. 5.90 × 10 −4 s −1 This study. The phosphorylated transcription factor can undergo de-phosphorylation such that This system of reactions leads to the following additional (for TF Atf3 , DNA Atf3 and Atf3·DNA Atf3 ) and revised (for Atf3 and TF Atf3 ·DNA Atf3 ) ODEs with at t = 0 and all other initial conditions are zero.
Conservation of D A (adding equations (42), (43), (45) and (46), integrating with respect to time and applying the initial conditions) leads to A similar relationship holds for F Assuming equation (46) is quasi-steady and utilising equation (48) leads to which on re-arranging gives where K 11 = k −11 /k 11 . This leads to the revised system of equations with the initial conditions

S2.1 Parameterisation, model solutions and sensitivity analysis
This model requires three new parameters: (i) the association constant for which Atf3 protein associates with Atf3 DNA (K 11 ); and (ii) the forward and reverse rates of the Atf3 transcription factor for Atf3 DNA. These were informed as follows. K 11 was assumed to have an initial value of K 11 = 0.1nM, a value in line with that previously found in (S6). This value was varied in order to understand how Atf3 association for its DNA affected the Atf3 mRNA and Atf3 concentration profiles as discussed in the main text. The ratio of λ −1 /λ 1 was assumed to have a similar value to K 11 , but a slow rate of reversal, such that λ −1 = 1×10 −4 /s was chosen. This leads to a value of λ 1 = 1×10 5 (M s) −1 . Variations to λ 1 and λ −1 (informed by a sensitivity analysis) showed that decreasing this ratio 2-fold had no effect on the results shown in Figure S3 in Text S1. Increasing the ratio 2-fold decreased the amount of bound Atf3 transcription factor and thus the amount of Atf3, which increased the amount of Egr1 which does agree with the experimental data shown in Figure 1B. The effect of increasing the rate of association of Atf3 for Atf3 DNA (K 11 ) is shown in Figures 1B and 1C of the main text.

S3 Model Extension 2 -RSK-P inhibits Atf3 DNA transcription
Here Pathway 2 in Figure S1 in Text S1 acts and leads to RSK being phosphorylated by ERK-P ERK-P + RSK k 9 → ERK-P + RSK-P.
Two transcription factors then compete for sites on the Atf3 DNA, that which binds to Atf3 and is phosphorylated by ERK-P (as per equations (38) and (39)), and RTF RTF + DNA Atf3 which is subsequently phosphorylated by RSK-P RTF · DNA Atf3 + RSK-P k 10 k −10 RTF-P · DNA Atf3 + RSK-P.
The additional and revised governing ODEs for these reactions are given by (70) such that and likewise that of RTF We also note that equation (49) still holds.
Bringing these results together, the revised systems of governing equations is given by with the initial conditions

S3.1 Parameterisation, model solutions and sensitvity analysis
This model extension leads to the new parameters R 0 (initial concentration of RSK), k 9 (rate of RSK phosphorylation by ERK-P), d 10 (rate of RSK-P de-phosphorylation), λ 2 and λ −2 (the rates of association and disassociation of RTF for Atf3 DNA) and k 10 and k −10 (the rates of phosphorylation and de-phosphorylation of RTF bound to Atf3 DNA). A fit-by-eye to the RSK data shown in Figure 3C yielded k 9 = 1 × 10 5 (Ms) −1 and d 10 = 5.9 × 10 −4 s −1 . We assume that λ 2 =λ 1 and λ −2 =λ −1 , k 10 = k 5 , k −10 = k −5 and R 0 = E 0 (the total amount of RSK is assumed equal to that of ERK). Solutions to the model equations are shown in Figure S4 in Text S1. The effect of varying the competition between TF and RTF for Atf3 DNA on Atf3 mRNA, Atf3 protein and Egr1 mRNA levels (sensitivity analysis) are discussed in the main text.

S4 Model extension 3 -ITF suppresses Atf3 transcription
In this case Pathway 3 in Figure S1 in Text S1 acts and leads to ITF protein being transcribed by RSK-P. This pathway begins by RSK-P phosphorylating the transcription factor bound ITF DNA such that which subsequently produces ITF mRNA or de-phosphorylates The ITF mRNA is subsequently translated to ITF protein whereby ITF protein subsequently binds to Atf3 DNA to inhibit transcription Finally the ITF mRNA and protein each degrade Applying the law of mass action to this series of reactions leads to As with the previous pathways we first observe that addition of equations (99), (101), (102) and (103) leads to conservation of Atf3 DNA such that Likewise conservation of the transcription factor bound to Atf3 DNA holds as defined by equation (49) and addition of equations (95) and (96) gives Assuming equation (103) is quasi-steady leads to Utilising these results leads to the simplified system of with the initial conditions

S4.1 Parameterisation, model solutions and sensitivity analysis
We require six further parameters for this model: (i) the disassociation constant of ITF protein for Atf3 DNA (K 15 ), which we assume is equal to that of RSK-P for Atf3 DNA (K 15 = K 10 ); (ii) ITF transcription (k 13 ) which we assume equal to that of Atf3 transcription (k 13 = k 6 ); (iii) translation rate of ITF protein (k 14 ) which we assume equal to that of Atf3 translation (k 14 = k 8 ); (iv) bound ITF transcription factor dephosphorylation (d p3 ), which we assume equal to that of Atf3 (d p3 =d p2 ); (v) ITF mRNA degradation (d 4 ) which we assume equal to Atf3 mRNA degradation (d 4 = d 2 ); and (vi) ITF degradation (d 5 ) which we assume equal to that of Atf3 (d 5 = d 3 ).
Solutions to the model equations are shown in Figure S5 in Text S1 and are discussed in the main text. A sensitivity analysis of the model (methodology detailed in Section S1.2.1) revealed that increasing or decreasing K 15 had no affect on the model solutions; the concentration in ITF is only briefly high enough to impeded Atf3 transcription (see Figure  S5(d) in Text S1) and varying the disassociation rate 2 orders of magnitude up or down was not enough to affect the overall Atf3 response. Increasing/increasing k 13 led to increases in ITF levels, but this do not affect the Atf3 mRNA or Atf3 profiles and subsequently the Egr1 mRNA profile remained the same. Varying the remaining parameters, k 14 , d 4 , d 5 and d p3 altered the amounts of ITF subsequently produced (mRNA and protein). Both of these variations were again not high to alter the Atf3 mRNA or Atf3 profiles and subsequently Egr1 mRNA given the transient nature of ITF binding to Atf3 DNA.

S5 Model Extension 4 -RSK-P driven microRNA expression increases Atf3 mRNA degradation and/or inhibits protein synthesis
Here we consider Pathway 4 in Figure S1 in Text S1 in which microRNA expression is facilitated by the phosphorylated RSK. The transcribed microRNA is then free to bind to Atf3 RNA whereby it forms an inactive complex, thus reducing the quantity of mRNA available for Atf3 translation.
The phosphorylated transcription factor subsequently produces premiRNA or de-phosphorylates miRNA is produced from the premiRNA premiRNA k 18 → miRNA (128) and both subsequently degrade Finally the miRNA can bind to the Atf3 RNA to form an inactive complex which inhibits Atf3 RNA translation miRNA + mRNA Atf3 and the complex can subsequently be degraded such that The additional ODEs describing these reactions are given by From this system of equations we see that the quantity of miDNA transcription factor is conserved (summation of equations (132) and (133)) This allows us to reduce the system of governing ODEs for Pathway 4 to with the initial conditions

S5.1 Parameterisation, model solutions and sensitivity analysis
This extension to the model introduces eight new parameters. A number of these are currently unknown experimentally and we hence proceed to inform them as follows.
k 16 , k −16 -We assume the rate of phosphorylation for the transcription factor bound miDNA is equal to that of the transcription factors bound to Egr1 and Atf3 by ERK-P (k 16 = k 3 ). The reverse rate of phosphorylation (k −16 ) is assumed equal to k −3 .
k 17 -We assume the rate of miDNA transcription is the same as that of Atf3 mRNA transcription, i.e. k 17 = k 6 .
k 18 -The rate at which premiRNA is converted to miRNA is assumed to be 0.5/s (double the rate of Atf3 translation).
k 19 -The rate at which miRNA associates with Atf3 mRNA is assumed to be 1×10 5 (Ms) −1 (equivalent to the rate we have set for Atf3 protein binding Egr1 DNA).
k −19 -The reverse rate at which miRNA associates with Atf3 mRNA is assumed to be 5×10 −5 /s -the value that was found to give the best qualitative fit to the experimental data.
d 6 , d 7 and d 8 -We assume the rate of decay of premiRNA (d 6 ) and miRNA (d 7 ) takes approximately 24 hours, which is line with recently published data (S5; 10). This leads to a decay rate of d 6 = 8.02×10 −6 = d 7 . The rate of decay of the complex miRNA·mRNA Atf3 is assumed the same as that of Atf3 mRNA, i.e. d 8 = d 3 .
d p4 -This is assumed equal to that of phosphorylated Atf3 bound transcription factor degradation (d p4 = d p2 ).
The effect of introducing the miRNA process on the expression levels of Atf3 mRNA is shown in Figure S6 in Text S1 and discussed in further detail in the main text. A sensitivity analysis of the model parameters (see Section S1.2.1 for methodology) revealed that when k 16 was increased/decreased the Atf3 mRNA, Atf3 and Egr1 profiles remained the same, but less/more Atf3 was produced. This is because increasing/decreasing RSK-P phosphosrylation of the transcription factor bound to the miDNA produces more/less miRNA which leads to less/more Atf3 mRNA and subsequently Atf3. However, this affect is not so large over the two orders of magnitude variation in k 16 so as to dramatically affect the concentration of Atf3 mRNA, Atf3 and Egr1 mRNA. The same effect was observed when altering k −16 . It was found that varying k 19 and k −19 had the same effect due to the immediate association between miRNA and Atf3 mRNA in increasing/decreasing levels of Atf3 mRNA.
Increasing k 17 decreased the amount of Atf3 mRNA and thus Atf3 which subseqently led to an increase in Egr1 mRNA, such that it no longer matched the experimental data. A similar result was found when k 17 was decreased meaning that the Atf3 mRNA profile no longer fitted the experimental data. In each case the Egr1 mRNA profile remained the same. Similar effects were observed when k 18 was subsequently increased or decreased. When increasing and decreasing the rates of degradation d 6 and d 7 the solutions for Atf3 mRNA, Atf3 and Egr1 mRNA remained the same; the concentration of Atf3 mRNA and Atf3 either rose or fell, but not significantly. Increasing or decreasing d 8 had enough effect to allow improved fitting of the Atf3 mRNA solution with the experimental data (see main text for discussion). Increasing d p4 had more of an effect on the model -data fitting as it more strongly influences how much miRNA is produced early on in the signalling cascade; in this case the Atf3 mRNA profile no longer matches the experimental data as not enough miRNA was produced to reduce the Atf3 mRNA levels. However, the Egr1 mRNA remained the same. Decreasing d p4 had no appreciable effect on the solutions.