First remove all existing global variables, and define the GK function

Remove["Global`*"]

G[u_,v_,J_,K_]:=(2u K)/(v-u+v J+u K+(Sqrt[(((v-u+v J+u K))^2)-4(v-u) u K]))

Model Figure 1F - note that k3 and k4 in Parameter set 1f have to be interchanged

Kinetic model of the network

rates1f={
v1->k0+k1 S,
v2->k2 R+k2p E R
}/.E->G[k3,k4 R,J3,J4]
param1f={
k0->0,
k1->1/20,
k2->1/10,
k2p->1/2,
k3->1/5,
k4->1,
J3->1/20,
J4->1/20}
odes1f={
R'[t]==Evaluate[(v1-v2)/.rates1f/.{R->R[t]}]
}

Plot[G[k3,k4 R,J3,J4]/.param1f,{R,0,1.5},PlotRange->All]

Manipulate[Plot[Evaluate[{v1,v2}/.rates1f/.param1f/.S->s],{R,0,1.5},AxesLabel->{"R","Rate (dR/dt)"},PlotRange->{0,.1}],{{s,1},0,2,Appearance->"Open"}]

Investigate how the number of intersections (stationary points) depends upon the value of s. What are the critical values of s where this number of intersections changes?

Now let us plot v1 - v2 as function of R and S, and see where it cuts the zero plane

Plot3D[Evaluate[{v1-v2,0}/.rates1f/.param1f],{R,0,1.5}, {S,0,2},AxesLabel->{R,S,R'},PlotStyle->{Red,Gray}]

Now let us calculate the intersection with the zero plane and calculate Rss as function of S

Sssrules=Solve[0==(v1-v2)/.rates1f/.param1f,S]//Flatten

Note the many complex numbers and the (□^(1/3))'s
Note that no Flatten was used, to arrive at a list with 3 solutions

Sss=S/.Sssrules

Plot[Sss,{R,0,1.5},PlotStyle->Thickness[0.01],AxesLabel->{"R","Sss"}]

ParametricPlot[{Sss,R},{R,0,1.5},PlotStyle->Thickness[0.01],AxesLabel->{"S","Rss"}]

plot1f[tend_]:=Module[{rule1f,rule1fback,plot1,plot2,plot3},rule1f=NDSolve[{odes1f/.param1f/.S->t,R[0]==0},{R},{t,0,tend}]//Flatten;Rtend=R[tend]/.rule1f;plot1=Plot[R[t]/.rule1f,{t,0,tend},AxesLabel->{"t","R"},PlotRange->All];rule1fback=NDSolve[{odes1f/.param1f/.S->(2tend-t),R[tend]==Rtend},{R},{t,tend,2tend}]//Flatten;plot2=Plot[Evaluate[R[t]/.rule1fback],{t,tend,2tend},PlotRange->All,AxesLabel->{"t","R"},PlotStyle->{Red,Dashed,Thick}];R2tend=R[2tend]/.rule1fback;rule2eback=NDSolve[{odes1f/.param1f/.S->0,R[2tend]==R2tend},{R},{t,2tend,5tend}]//Flatten;plot3=Plot[Evaluate[R[t]/.rule2eback],{t,2tend,5tend},PlotRange->All,AxesLabel->{"t","R"},PlotStyle->{Green,Dashed,Thick}];{Show[plot1,plot2,plot3,PlotRange->All],Plot[Piecewise[{{t,0<t<tend},{(2tend-t),tend<t<2tend},{0,t>2tend}}],{t,0,5tend},PlotRange->All,AxesLabel->{"t","S"}]}]

Manipulate[plot1f[tend],{tend,0.001,10}]

Now determine the critical S values !

plot1f2[tend_,s_]:=Module[{rule1f,rule1fback,plot1,plot2},rule1f=NDSolve[{odes1f/.param1f/.S->s,R[0]==0},{R},{t,0,tend}]//Flatten;Rtend=R[tend]/.rule1f;plot1=Plot[R[t]/.rule1f,{t,0,tend},AxesLabel->{"t","R"},PlotRange->All];rule1fback=NDSolve[{odes1f/.param1f/.S->(s+tend-t),R[tend]==Rtend},{R},{t,tend,s+tend}]//Flatten;plot2=Plot[Evaluate[R[t]/.rule1fback],{t,tend,s+tend},PlotRange->All,AxesLabel->{"t","R"},PlotStyle->{Red,Dashed,Thick}];R2tend=R[s+tend]/.rule1fback;rule2eback=NDSolve[{odes1f/.param1f/.S->0,R[s+tend]==R2tend},{R},{t,s+tend,s+2tend}]//Flatten;plot3=Plot[Evaluate[R[t]/.rule2eback],{t,s+tend,s+2tend},PlotRange->All,AxesLabel->{"t","R"},PlotStyle->{Green,Dashed,Thick}];Show[plot1,plot2,plot3,PlotRange->All]]

Manipulate[plot1f2[tend,s],{{tend,20},0.001,100},{s,0.001,2}]

To determine stability we must linearize around a stationary point

FRS=(v1-v2)/.rates1f/.param1f

Rsspoint[s_]:=R/.Solve[0==FRS/.S->s,R]

Rsspoint[1.5]

D[FRS,R]/.R->Rsspoint[1.5][[1]]

D[FRS,R]/.R->Rsspoint[1.5]

negative = stabile, positive = instabile

to find the linearized equation of FRS we can use Series and Normal

nFRS=Normal[Series[FRS,{R,Rsspoint[1.5][[1]],1}]]/.S->1.5//Factor

control that this is indeed still a stationary point

nFRS/.R->Rsspoint[1.5][[1]]

and this function summarizes the linearized equation (s) for arbitrary S

linearize[FRS_,s_]:=Module[{Rsspoints},Rsspoints=R/.Solve[0==FRS/.S->s,R];Map[Normal[Series[FRS,{R,#,1}]]&,Rsspoints]/.S->s//Factor]

linearize[FRS,1.5]

Model Figure 1G

Kinetic model of the network

rates1g={
v1->k0 E,
v2->k2 S R
}/.E->G[k3,k4 R,J3,J4]
param1g={
k0->1,
k2->1,
k3->1/2,
k4->1,
J3->1/100,
J4->1/100};
odes1g={
R'[t]==Evaluate[(v1-v2)/.rates1g/.{R->R[t]}]
}

Manipulate[Plot[Evaluate[{v1,v2}/.rates1g/.param1g/.S->s],{R,0,1},AxesLabel->{"R","Rate (dR/dt)"},PlotRange->{0,1.2}],{{s,1},0,2}]

Now let us plot v1 - v2 as function of R and S, and see where it cuts the zero plane

Plot3D[Evaluate[{v1-v2,0}/.rates1g/.param1g],{R,0,1}, {S,0,2},AxesLabel->{R,S,R'},PlotStyle->{Red,Gray}]

Now let us calculate the intersection with the zero plane and calculate Rss as function of S

Rssrules=Solve[0==(v1-v2)/.rates1g/.param1g,R]

Plot[Evaluate[R/.Rssrules/.param1g],{S,0,2},PlotStyle->Thickness[0.01],AxesLabel->{"S","Rss"}]

The middle solution turns out to be the correct solution. Most probably the reason for the spurious solutions is related to the fact that Mathematica treats all variables as complex, and the (□^(1/3))'s.

Exercise Model Figure 1E, to be handed in

First define the Kinetic model of the network

now plot v1 and v2, reproducing the third column of Fig .1 e

Now plot v1 - v2 as function of R and S, and see where it cuts the zero plane

Now calculate the intersection with the zero plane and calculate Rss as function of S

And plot this result

For the dynamics we choose a simple signal S = t, for t from 0 till 40, using R(0)=0.

Write a caption for this plot

Make a plot of R (t), for t from 0 till 40

Calculate the value of R (t) at t = 40, and call the object R40

and now let the signal decrease back to zero again S=80-t , demonstrating the irreversible switch. calculate for t from 40 till 80, using R(40)=R40

Make a plot of R (t), for t from 40 till 80

Calculate the value of R (t) at t = 80, and call the object R80

and now let the signal be zero S=0, calculate for t from 80 till 120, using R(80)=R80

Make a plot of R (t), for t from 80 till 120

Make a plot of the signal S for t from 0 till 120, and show all plots.

Write a caption for this plot

Now copy the above function plot1f to a new function plot1fe and modify it such that it  automates the previous process

Now determine the critical S value ! What is it?

And summarize the linearized equation (s) for S values between 0 and 15 and comment on the outcome for the different regimes. What is the stability of the stationary points?