the nanomolar range and the level of H
2
O
2
is approximately 50-fold higher than the
concentration of superoxide. Finally, we use Plot to produce a graphic showing the
time course of the variable concentrations. We specify several options such as axes la-
bels and colors to make the plot more informative.
sol = NDSolve[{eq1, eq2, O2[0]==0, H2O2[0]==0}/.par, {O2, H2O2}, {t,0,0.01}]
{O2[0.01], H2O2[0.01]}/.sol ⇒ {{4.125*10
–11
, 1.84552*10
–9
}}
Plot[Evaluate[{O2[t], H2O2[t]/50}/.sol], {t,0,0.01}, PlotRange->All,
PlotStyle->{Hue[0.9],Hue[0.6]}, AxesLabel->{“time”,“concentration”}];
14.1.1.2 Matlab Example
In Matlab the routine that solves the ODE systems requires as one of its arguments
a function that evaluates the right-hand side of the ODE system for given times and
variable values. Because the example system is so small, we can define this function
as an inline function and avoid writing a separate M-file. Next the options for the
ODE solver are defined. The absolute values of the solution are very small, and there-
fore the absolute tolerance has to be adjusted accordingly.
dydt = inline(‘[c1-c2*SOD*y(1);c2*SOD*y(1)-c3*cat*y(2)]‘,’t‘,’y‘,’tmp‘,’
SOD‘,’cat‘,’c1‘,’c2‘,’c3’);
options = odeset(‘OutputFcn’,[ ], ‘AbsTol’, 1e-30, ‘RelTol’, 1e-6);
Now ode45 is called, which solves the ODE system for the specified time span. In
addition to dydt, which evaluates the derivatives,we also supply the starting concentra-
tions of the variables and the numerical values for the constants that were used in
dydt. ode45 is only one of a whole set (ode45, ode23, ode113, ode15s, ode23s, ode23t,
and ode23tb) of possible solvers with different properties (differing in accuracy or suit-
ability for stiff problems). The function returns a vector, t, holding time points for
which a solution is returned and a matrix, y, that contains the variable values at these
time points. To find out which concentration exists at a given time, the function “in-
terp1” can be used, which also interpolates between existing time points if necessary.
Finally, the resulting time course is plotted and the axes are labeled appropriately.
[t, y] = ode45(dydt, [0,0.01], [0;0], options, 1e-5, 1e-5, 6.6e-7, 1.6e9, 3.4e7);
interp1(t,y,0.01) ⇒ 4.125*10
–11
, 1.84552*10
–9
plot(t, y(:,1), t, y(:,2)/50)
xlabel(‘time’)
ylabel(‘concentration’)
14.1.2
Gepasi
Matlab and Mathematica are huge and expensive general-purpose tools for mathe-
matical modeling. They can be used to model anything that can be modeled, but at
the cost of a steep learning curve. The opposite approach is used by specialized tools
422
14 Modeling Tools