Dynamic behavior of autocatalytic reaction networks
Authors' Data: Housam Binous & Nasri Zakia
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Reference :
E. Di Cera, P. E. Phillipson and J. Wyman, Limit-cycle oscillations and chaos in reaction networks subject to conservation of mass, Proc. Natl. Acad. Sci., Vol. 86, pp. 142-146, January 1989.

Period 1 behavior for m=0.5

sol[m_]:=NDSolve[{D[x1[t],{t,1}]== x1[t]((Subscript[K, 4]) x4[t]-(Subscript[K, 1])(m-x2[t])x2[t]),
D[x2[t],{t,1}]== x2[t]((Subscript[K, 1]) (m-x2[t])x1[t]-(Subscript[K, 2])x3[t]),
D[x3[t],{t,1}]== x3[t]((Subscript[K, 2]) x2[t]-(Subscript[K, 3])x4[t]),
D[x4[t],{t,1}]== x4[t]((Subscript[K, 3]) x3[t]-(Subscript[K, 4])x1[t]),x1[0]==0.97,x2[0]==0.01,x3[0]==0.01,x4[0]== 0.01},{x1[t],x2[t],x3[t],x4[t]},{t,0,10000},MaxSteps-> 100000]

(Subscript[K, 1])=10;(Subscript[K, 2])=1;(Subscript[K, 3])=1;(Subscript[K, 4])=1;

f=sol[0.5]

Plot[First[Evaluate[x3[t]/.f]],{t,5000,6000},PlotStyle-> RGBColor[1,0,0],PlotRange->All]

ParametricPlot[{First[Evaluate[x3[t]/.f]],First[Evaluate[x4[t]/.f]]},{t,5000,6000},PlotStyle-> RGBColor[1,0,0],PlotRange->{{0,1},{0,1}}]

Period 2 behavior for m=0.56

f=sol[0.56]

Plot[First[Evaluate[x3[t]/.f]],{t,5000,6000},PlotRange-> All,PlotStyle-> RGBColor[1,0,0]]

ParametricPlot[{First[Evaluate[x3[t]/.f]],First[Evaluate[x4[t]/.f]]},{t,5000,6000},PlotStyle-> RGBColor[1,0,0],PlotRange->{{0,1},{0,1}}]

Period 4 behavior for m=0.57

f=sol[0.57]

Plot[First[Evaluate[x3[t]/.f]],{t,5000,6000},PlotRange-> All,PlotStyle-> RGBColor[1,0,0]]

ParametricPlot[{First[Evaluate[x3[t]/.f]],First[Evaluate[x4[t]/.f]]},{t,5000,6000},PlotStyle-> RGBColor[1,0,0],PlotRange->{{0,1},{0,1}}]

Period 8 behavior for m=0.572

f=sol[0.572]

Plot[First[Evaluate[x3[t]/.f]],{t,5000,6000},PlotStyle-> RGBColor[1,0,0],PlotRange-> All]

ParametricPlot[{First[Evaluate[x3[t]/.f]],First[Evaluate[x4[t]/.f]]},{t,5000,6000},PlotStyle-> RGBColor[1,0,0],PlotRange->{{0,1},{0,1}}]

Chaotic behavior for m=0.575

f=sol[0.575]

Plot[First[Evaluate[x3[t]/.f]],{t,5000,6000},PlotRange-> All,PlotStyle-> RGBColor[1,0,0]]

ParametricPlot[{First[Evaluate[x3[t]/.f]],First[Evaluate[x4[t]/.f]]},{t,5000,6000},PlotStyle-> RGBColor[1,0,0],PlotRange->{{0,1},{0,1}}]

Bifurcation diagram

p=0; μ=0.5;

While[p<1001  && μ< 0.578 ,{μ=0.5+0.0002*p,f=sol[μ],If[μ<0.56,p+=10,p+=1],Print[μ," ",p],i=0,k=0,While[i<2001,{t0=5000+i,i+=10,solt[i]=Quiet[FindRoot[D[First[x3[t]/.f],{t,1}]==0,{t,t0}]],If[5000<(t/.solt[i])<7000 &&Quiet[(D[First[x3[t]/.f],{t,2}]/.solt[i])<0],{pt[k]={μ,First[x3[t]/.f]/.solt[i]},k+=1}]}],tbl[p]=Table[pt[i],{i,1,k-1}]}]

ListPlot[Join[Table[tbl[i],{i,10,300,10}],Table[tbl[i],{i,300,380,1}]],PlotRange-> All,AspectRatio->1,Frame->True,PlotStyle-> Table[{PointSize[0.0059],RGBColor[1,0,0]},{i,1,101}]]

ListPlot[Table[tbl[i],{i,301,380,1}],AspectRatio->1,AxesOrigin->{0.57,0.4},Frame-> True,PlotStyle-> Table[{PointSize[0.0059],RGBColor[1,0,0]},{i,1,1001}]]