<< myPadRight.m;
<< 
        cellerator.m;
<< 
        kMech.m;
Off[General::"spell1"];
Off[General::"spell"];

(* 
      Start Transcription/
          RNA Splcing Coupling Simulation *) 

 (* NE1 :  RNA polymerase;  NEs :  spliceosome;  
      NEd :  degradation enzymes *) 

mRNASplcing = 
      Union[
Enz[(Overscript[({DNA, NTP} ⇄ {premRNA, DNA}), NE1]), 
          BiBi[kf1, kr2, 
            kcat3]], 
{{premRNA + 
                NEs ⇄ $Complex$NEs$premRNA$, kf4, 
            kr5}}, 
{{$Complex$NEs$premRNA$ → NEs + mRNA, 
            kcat6}}, 
{{(Overscript[(premRNA ⇄ NMP),            NEd]), kf7, kr8, 
            kcat9}}, 
{{(Overscript[($Complex$NEs$premRNA$ ⇄ NMP),                                                      NEd]), 
            kf10, kr11, 
            kcat12}}, 
{{(Overscript[(mRNA ⇄ NMP), NEd]), kf13, kr14, kcat15}}
]

{myODEs,myVars}=interpret[mRNASplcing]

Lamda  =  1;

    
\(myKs = {

KmNE1$DNA = 40; 
        KmNE1$NTP = 1000; 
        kcat3 = 0.1; 
kf1 -> 
          Kf2S[KmNE1$DNA, KmNE1$NTP, kcat3, 
            Lamda], 
kr2 -> 
          Kr[kcat3, 
            Lamda], 

KmNEs$premRNA = 
          3.5;  (*\(3.5;\)*) kcat6 = 0.03; 
kf4 -> 
          Kf[KmNEs$premRNA, kcat6, Lamda], 
kr5 -> 
          Kr[kcat6, Lamda], 

\ \ \ \ KmNEd$premRNA = 
          2000000; 
        kcat9 = 30000; 
kf7 -> 
          Kf[KmNEd$premRNA, kcat9, Lamda], 
kr8 -> 
          Kr[kcat9, 
            Lamda], 

KmNEd$Complex$NEs$premRNA = 30000000; 
        kcat12 = 30000; 
kf10 -> 
          Kf[KmNEd$Complex$NEs$premRNA, kcat12, 
            Lamda], 
kr11 -> 
          Kr[kcat12, 
            Lamda], 

KmNEd$mRNA = 
          30000000; \ 
        kcat15 = 30000; 
kf13 -> 
          Kf[KmNEd$mRNA, kcat15, Lamda], 
kr14 -> 
          Kr[kcat15, Lamda]
}\)

myICs = {
DNA[0] == 
            12, 
NTP[0] == 
            25000, 
premRNA[0] == 
            0, 
mRNA[0] == 
            0, 
NMP[0] == 
            0, 

NE1[0] == 
            0.025, 
NEs[0] == 
            5, 
NEd[0] == 
            3.5, 

$Complex$NE1$DNA$NTP$[
              0] == 
            0, 
$Complex$NEs$premRNA$[0] == 
            0, 
$Complex$premRNA$NEd$[0] == 
            0, 
$Complex$$Complex$NEs$premRNA$$NEd$[
              0] == 
            0, 
$Complex$mRNA$NEd$[0] == 
            0
};

tmax = 165;
mySolution = 
      NDSolve[Join[myODEs /. myKs, myICs], myVars, {t, 0, tmax}, 
        MaxSteps -> 3000]

orderedVars={DNA,NTP,premRNA,mRNA,NMP,NE1,NEs,NEd,$Complex$NE1$DNA$NTP$,
$Complex$NEs$premRNA$,$Complex$premRNA$NEd$,$Complex$mRNA$NEd$};

Show[GraphicsArray[
    Partition[
      Map[Plot[#[t]/.mySolution,{t,0,tmax},PlotLabel->#,
            PlotRange->All,DisplayFunction->Identity]&,orderedVars],3]
    ]
  ];

FormalPlot[varName_, displayTime_, 
          YplotRange_] := 
Plot[
          varName[t] /. mySolution, {t, 0, displayTime}, 
          PlotRange -> YplotRange, Background -> GrayLevel[1], 
          TextStyle -> {FontWeight -> "Bold", 
              FontSize -> 12, 
              FontFamily ->  "Times New Roman"}, 
          PlotStyle -> {Thickness[0.008]}, 
          AxesStyle -> {Thickness[
                0.0006]}];
DashPlot[varName_, displayTime_, 
          YplotRange_] := 
Plot[
          varName[t] /. mySolution, {t, 0, displayTime}, 
          PlotRange -> YplotRange, Background -> GrayLevel[1], 
          TextStyle -> {FontWeight -> "Bold", 
              FontSize -> 12, 
              FontFamily ->  "Times New Roman"}, 
          PlotStyle -> {Thickness[0.008], Dashing[{0.02, 0.025}]}, 
          AxesStyle -> {Thickness[
                0.0006]}];
DotPlot[varName_, displayTime_, 
          YplotRange_] := 
Plot[
          varName[t] /. mySolution, {t, 0, displayTime}, 
          PlotRange -> YplotRange, Background -> GrayLevel[1], 
          TextStyle -> {FontWeight -> "Bold", 
              FontSize -> 12, 
              FontFamily ->  "Times New Roman"}, 
          PlotStyle -> {Thickness[0.008], 
              Dashing[{0.001, 0.001, 0.001, 0.015}]}, 
          AxesStyle -> {Thickness[0.0006]}];

totalPremRNA[t_] = 
        premRNA[t] + $Complex$NEs$premRNA$[t] + $Complex$premRNA$NEd$[
            t];

P1 = DotPlot[mRNA, tmax, {0, 0.1}];
P2 = 
        Plot[totalPremRNA[t] /. mySolution, {t, 0, tmax}, 
          PlotRange -> All, 
          DisplayFunction -> Identity];
Show[P1, P2, PlotRange ->  {0, 0.1}, 
        DisplayFunction -> $DisplayFunction];

<< Graphics`MultipleListPlot`

  

  
mRNA = {{15, 0.0019}, {40, 0.0081}, {65, 0.024}, {90, 0.041}, {115, 
            0.057}, {140, 0.0725}, {165, 
            0.084}};

  
premRNA = {{15, 0.0139}, {40, 0.04}, {65, 0.056}, {90, 0.053}, {115, 
            0.052}, {140, 0.055}, {165, 
            0.056}};

P3 = 
        MultipleListPlot[mRNA, premRNA, 
          SymbolShape -> {PlotSymbol[Triangle, 5], 
              PlotSymbol[Box, 2.5]}, Background -> GrayLevel[1], 
          TextStyle -> {FontWeight -> "Bold", 
              FontSize -> 12, 
              FontFamily ->  "Times New Roman"}, 
          PlotStyle -> {Thickness[0.008], Dashing[{0.02, 0.025}]}, 
          AxesStyle -> {Thickness[0.0006]}];

Show[P3, P1, P2, PlotRange ->  {0, 0.1}, 
        DisplayFunction -> $DisplayFunction];