The Steinmetz-Larter example

In the peroxidase-oxidase reaction a peroxidase catalyzes an aerobic oxidation:

\begin{displaymath}
2YH_2+O_2+2H^+\rightarrow 2YH^+ + 2H_2O,
\end{displaymath}

where $YH_2$ ($NADH$) is a general electron donor. Under the proper conditions this reaction yields oscillations in the concentrations of $YH_2$, $O_2$ and various reaction intermediates. Steinmetz and Larter [Steinmetz and Larter 1991] derived the following system of four coupled nonlinear differential rate equations:
\begin{displaymath}
\left\{\begin{array}{rcl}
\dot A&=& -k_1 ABX- k_3 ABY+ k_7...
...,\\
\dot Y&=& -k_3 ABY+2k_2X^2-k_5 Y,
\end{array}\right.
\end{displaymath} (32)

where $A,B,X,Y$ are state variables and $k_1$, $k_2$, $k_3$, $k_4$, $k_5$, $k_6$, $k_7$, $k_8$, and $k_{-7}$ are parameters. Although highly simplified, the model is able to reproduce the three modes of simple, chaotic, and bursting oscillations found in experiments. State and parameter values of a point on a NS cycle in (32) are given in Table 4. The normal form coefficient of the NS cycle is $-1.406017e-006$. Since it is negative, we expect stable tori nearby.


Table 4: State and parameter values of a point on an NS- cycle in the Steinmetz-Larter model.
Variable Value Parameter Value Parameter Value
$A$ 1.8609653 $k_1$ 0.1631021 $k_5$ 1.104
$B$ 25.678306 $k_2$ 1250 $k_6$ 0.001
$X$ 0.010838258 $k_3$ 0.046875 $k_7$ 0.71643356
$Y$ 0.094707061 $k_4$ 20 $k_{-7}$ 0.1175
$k_8$ 0.5


To start a time integration from this NS point, with a slightly supercritical parameter value, namely $k_7=0.7167$, we can use the commands

X0=[1.8609653;25.678306;0.010838258;0.094707061];
P0=[0.1631021;1250;0.046875;20;1.104;0.001;0.71643356;0.1175;0.5];
P0(7)=0.7167;
hls=feval(@SteLar);
options=odeset('RelTol',1e-8);
[t,y]=ode45(hls{2},[0,3000],X0,options,0.1631021,1250,0.046875,20,1.104,
0.001,0.7167,0.1175,0.5);
size(t)
size(y)
plot(t(100000:339000),y(100000:339000,2));
'press any key'
pause
plot(y(100000:339000,3),y(100000:339000,2));

After a transient the time-series exhibits modulated oscillations with two frequencies near the original limit cycle (see Fig. 7). This is a motion on a stable two-dimensional torus that arises from the Neimark-Sacker bifurcation. The commands can be executed by running the testrun testtimeinteg.m. The total number of computed points is 339161.

Figure 7: Movement on a stable torus in a 4-dimensional space as a function of time in the $B-$ direction and as projected on the $(B,X)$ phase plane.
\begin{figure}
\psfig{figure=ex/moduloscill.eps,height=4cm}\hfill \psfig{figure...
...a) Modulated oscillations.\hfill(b) Orbits on a stable 2-torus.}
\end{figure}

The testrun testtimeintegJacobian.m

X0=[1.8609653;25.678306;0.010838258;0.094707061];
P0=[0.1631021;1250;0.046875;20;1.104;0.001;0.71643356;0.1175;0.5];
P0(7)=0.7167;
hls=feval(@SteLar);
options=odeset('RelTol',1e-8);
options = odeset('Jacobian',hls(3),'JacobianP',hls(4),
'Hessians',hls(5),'HessiansP',hls(6));
[t,y]=ode45(hls{2},[0,3000],X0,options,0.1631021,1250,0.046875,20,
1.104,0.001,0.7167,0.1175,0.5);
size(t)
size(y)
plot(t(100000:325000),y(100000:325000,2));
'press any key'
pause
plot(y(100000:325000,3),y(100000:325000,2));

is a small modification of testtimeinteg.m in which the solver has access to the Jacobian matrix. The graphical output is similar, but the number of computed points is now 328197.