ABSTRACT (original article): In the present paper, we present the solution of four problems drawn from the chemical and biochemical engineering field of study. These problems illustrate various important aspects of nonlinear dynamics such as limit cycles, quasi-periodic and chaotic behaviors, time series and phase portraits, power spectra, the time-delay reconstruction diagrams, Hopf bifurcation, bifurcation diagrams, steady state multiplicity… Among a large number of examples, we have selected the following case studies because they illustrate the basic nonlinear dynamics concepts listed previously: (i) The glycolytic oscillator model, first suggested by Sel’kov in 1968 in order to elucidate the mechanism that living cells use to obtain energy from sugar breakdown; (ii) The Oregonator model derived by Field, Körös, and Noyes in the early 1970s, which elucidates the famous oscillatory behavior observed in the Belousov–Zhabotinsky reaction; (iii) The steady state multiplicity in a biochemical reactor for both the Monod and substrate inhibition kinetics; and (iv) The three-variable autocatalator, first proposed by Peng, Scott, and Showalter in 1990. Finally, we conclude this paper by sharing our experience teaching this subject to graduate students at King Fahd University of Petroleum & Minerals (KFUPM). CITATION (original article): Binous H., Bellagi A., Introducing nonlinear dynamics to chemical and biochemical engineering graduate students using MATHEMATICA©, Computer Applications in Engineering Education, 2019; 27: 217–235. https://doi.org/10.1002/cae.22070
In the study of nonlinear dynamics, it is useful to first introduce students to simple systems that exhibit periodic behavior as a consequence of a Hopf bifurcation. For instance, the glycolytic oscillator model has this feature. This model was first suggested by Sel’kov in 1968 in order to elucidate the mechanism that living cells use to obtain energy from sugar breakdown. This two-dimensional nonlinear and autonomous system, is given by
​
​

x
=
f
1
(x,y)=-x+ay+
2
x
y
​

​

y
=
f
2
(x,y)=b-ay-
2
x
y
These equations describe the autocatalytic reaction of two intermediate species
x
and
y
in an isothermal batch reactor, when the system is far from equilibrium [1]. In this context, the steady state referred to below is a pseudo steady state, and is applicable when the precursor reactant is slowly varying with time.
The unique steady state is given by
x
S
=b
and
y
S
=b/(a+
2
b
)
. This steady state is at the position of the green dot in the phase portrait diagram. It appears as the intersection of the dotted blue and green curves, which are the level curves given by
f
1
(x,y)=0
and
f
2
(x,y)=0
. The stability of steady state to small disturbances can be assessed by determining the eigenvalues of the Jacobian
J=
∂
f
1
/∂x
∂
f
1
/∂y
∂
f
2
/∂x
∂
f
2
/∂y
. A Hopf bifurcation occurs when a complex conjugate pair of eigenvalues crosses the imaginary axis.
We show the phase portrait with the vector field of directions around the critical point
(
x
S
,
y
S
)
. In addition, the eigenvalues of
J
, the trace
tr(J)
, the determinant
J
, and
2
Δ=tr(J))
-4J
are computed. The time series are also plotted in a separate diagram (the blue and red curves correspond to
x(t)
and
y(t)
, respectively. Finally, the real and imaginary parts of the two eigenvalues
λ
1
and
λ
2
are plotted (in blue and red, respectively) versus
b
for several values of the parameter
a
.
To observe a Hopf bifurcation, you can set
a=0.08
and vary
b
. As shown in this notebook, there are two Hopf bifurcation points. Indeed, for
b=0.3
, the steady state is a stable focus; for
b=0.5
, the trajectory is attracted to a stable periodic solution (called a limit cycle); and finally for
b=1.0
, the steady state is again a stable focus.
​
For
a=0.08
, the two Hopf bifurcation points are obtained at
b=0.34641
and
b=0.848528
.
Finally, for the last calculation (obtained for larger values of
a
, for instance,
a=0.65
), the Hopf bifurcation disappears and there are only transitions from a stable node to a stable focus and then back to stable node.
Two-dimensional nonlinear and autonomous system that describes the mechanism of glycolysis:
f1[x_,y_]:=-x+a*y+x^2y;​​f2[x_,y_]:=b-a*y-x^2*y;
Steady-states:
solSS=Flatten@NSolve[{-x+a*y+x^2y0,b-a*y-x^2*y0},{x,y}];
Jacobian matrix:
m1=Outer[D,{f1[x,y],f2[x,y]},{x,y}];
m=m1/.solSS;
Eigenvalues of Jacobian matrix:
Eig=Chop@N@Eigenvalues[m];
Initial condition:
x0={1.25,1.25};
a and b values:
a=0.08;b=0.30;
Solution of the dynamical system:
sol=NDSolve[{x'[t]-x[t]+a*y[t]+x[t]^2y[t],y'[t]b-a*y[t]-x[t]^2*y[t],x[0]x0[[1]],y[0]x0[[2]]},{x,y},{t,0,500},Method"StiffnessSwitching"];
Time series plot indicate damped oscillations - the eigenvalues are complex conjugates with negative real parts:
p1=Plot[Evaluate[{x[t],y[t]}/.sol],{t,0,80},PlotRangeAll,PlotStyle{Darker@Red,Thickness[0.007]},PlotPoints1000,AspectRatio1,FrameTrue,GridLinesAutomatic,FrameLabel{Style["time",16,Italic],Style["",16,Italic]},PlotPoints100,PlotLegends{Style["x",16,Italic],Style["y",16,Italic]}]
Out[]=
x
y
Eigenvalues:
zz=Chop@N@Eigenvalues[m];
Parametric plot:
p1=ParametricPlot[Evaluate[{x[t],y[t]}/.sol],{t,0,50},PlotRangeAll,PlotStyle{Darker@Red,Dashed,Thickness[0.007]},PlotPoints1000];​​xx=First[Cases[p1,x_LineFirst@x,Infinity]];
In[]:=
Show​​ContourPlot{-x+a*y+x^2y==0,b-a*y-x^2*y0},{x,0,3},{y,0,3},​​PlotLabelGridTranspose@​​"eigenvalues",IfIm[First@zz]0,NumberForm[zz,{4,2}],RowStyle["z",Italic]," = ",NumberForm[First@zz,{4,2}],", ",
Style["z",Italic]
,​​{"trace",NumberForm[Chop@N@Tr[m],{4,2}]},​​{"determinant",NumberForm[Chop@N@Det[m],{4,2}]},​​{"Δ",NumberForm[Chop@N@Tr[m]^2-4Det[m],{4,2}]},ItemSize({#,Automatic}&/@{10,4,4,4}),​​FrameLabel{Style["x",16,Italic],Style["y",16,Italic]},PlotPoints100,ContourStyle{{Green,Thick,Dashed,Opacity[0.99]},{Blue,Dashed,Thick,Opacity[0.99]}},Epilog{Green,PointSize[0.025],Point[{x,y}/.solSS],Cyan,Point[x0]},​​Graphics[{Red,Arrowheads[Table[0.06,{6}]],Dashed,Thickness[0.005],Arrow[Line[Partition[xx,Floor[Length[xx/10]]][[;;;;100]]]]}],​​StreamPlot[{-x+a*y+x^2y,b-a*y-x^2*y},{x,0,3},{y,0,3},StreamStyleDirective[Opacity[0.5]],StreamPoints50],​​ImageSize400{1,1},ImagePadding{{50,Automatic},{Automatic,Automatic}}
Out[]=
Let us determine the critical values of b when a=0.08
​
-- at these values both eigenvalues are purely imaginary shown below by the orange dots --
​
cyan dots are for b=0.3
b1c=0.3464101615137754;b2c=0.848528137423857;
Two-dimensional nonlinear and autonomous system that describes the mechanism of glycolysis:
f1b1[x_,y_]:=-x+a*y+x^2y;​​f2b1[x_,y_]:=b1-a*y-x^2*y;
Steady-states:
solSSb1=Flatten@Quiet@Solve[{-x+a*y+x^2y0,b1-a*y-x^2*y0},{x,y}];
Jacobian matrix:
m2=Outer[D,{f1b1[x,y],f2b1[x,y]},{x,y}]/.solSSb1;
Eigenvalues of Jacobian matrix:
{λ1,λ2}=Eigenvalues[m2];
Show[Plot[{Re[λ1],Re[λ2]},{b1,0,3},PlotStyle{{Red,Thick},{Blue,Dashed,Thick}},FrameTrue,FrameLabel{Style["b",Italic,16],Style["Re(λ)",16]},PlotRange{{0,3},{-5,1}},GridLinesAutomatic,AspectRatio1],​​Graphics[{Darker@Cyan,PointSize[0.03],Point[{b,Re[Eig[[1]]]}],Point[{b,Re[Eig[[2]]]}],Orange,PointSize[0.032],Point[{b1c,0}],Point[{b2c,0}]},ImagePadding{{40,Automatic},{40,Automatic}}],​​ImageSize.9×400{1,1}]
Show[​​Plot[{Im[λ1],Im[λ2]},{b1,0,3},PlotStyle{{Red,Thick},{Blue,Dashed,Thick}},FrameTrue,FrameLabel{Style["b",Italic,16],Style["Im(λ)",16]},PlotRange{{0,3},{-3,3}},GridLinesAutomatic,AspectRatio1],​​Graphics[{Darker@Cyan,PointSize[0.03],Point[{b,Im[Eig[[1]]]}],Point[{b,Im[Eig[[2]]]}]},ImagePadding{{40,Automatic},{40,Automatic}}],​​ImageSize.9×400{1,1}]
Out[]=
Determination of the critical values of b when a = 0.08
Color function for the contour plot given below:
Two-dimensional nonlinear and autonomous system that describes the mechanism of glycolysis:
Steady states:
Jacobian matrix:
Critical values of b when a=0.08:
We consider only positive b:
In[]:=
(*b=1.0 and a=0.08*)
In[]:=
(*b=0.4 and a=0.08*)
In[]:=
(*a=0.65 and b=0.05;*)

Acknowledgement

Some of the coding was written in collaboration with Prof. B. G. Higgins.

CITE THIS NOTEBOOK

The glycolytic oscillator mechanism​
by Housam Binous and Ahmed Bellagi
Wolfram Community, STAFF PICKS, April 8, 2024
​https://community.wolfram.com/groups/-/m/t/3154806