Finding a kinetic rate law using experimental 
data from a batch enzymatic reactor
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Acknowledgement :

This example was given by Kenneth J. Beers in Numerical Methods for Chemical Engineering, Applications in Matlab, Cambridge University Press, 2007.

Measured substrate concentration in batch enzymatic reactor

time = {30, 60, 90, 120, 180, 240, 300, 315, 330};

sub = {1.87, 1.73, 1.58, 1.43, 1.07, 0.63, 0.12, 0.04, 
          0.01};

TableForm[
      Transpose[{Join[{"time"}, time], 
          Join[{"Substarte"}, sub]}]]

mole balance on the substrate in a bioreactor containg the enzyme

myode[Vm_, Km_, Ksi_] := 
      NDSolve[{D[
              S[t], {t, 
                1}] ==  -mE/(α VR) (Vm S[
                      t])/(Km + S[t] + Ksi^-1 S[t]^2), 
          S[0] ==  2}, {S[t]}, {t, 0, 350}]

α = 10^6; VR = 100/1000; mE = 10;

Function that will be minimized (sum of squared differences between exprimental and predicted substrate concentration). We find Vm=211.184, Km=0.2258 and Ksi=0.5229.

sumsquare[Vm_?NumericQ, Km_?NumericQ, Ksi_?NumericQ] := 
      Sqrt[Sum[(sub[[i]] - 
                First[S[t] /. myode[Vm, Km, Ksi] /. 
                    t ->  time[[i]]])^2, {i, 1, 9}]]

FindMinimum[
      sumsquare[Vm, Km, Ksi], {Vm, 180, 220}, {Km, 0.1, 0.3}, {Ksi, 0.4, 
        0.6}]

Plotting predicted and experimental values of the substrate concentration versus time. We get perfect agreement.

plt1 = 
      Plot[Evaluate[S[t] /. myode[211.184, 0.225823, 0.522951]], {t, 0, 
          350}]

plt2 = 
      ListPlot[Table[{time[[i]], sub[[i]]}, {i, 1, 9}], 
        PlotStyle ->  {PointSize[0.025], RGBColor[1, 0, 0]}]

Show[plt1, plt2]

Plotting predicted and experimental values of the substrate concentration versus time. We get perfect agreement for the values found in the Book by Kenneth Beers using Matlab (Vm=200, Km=0.201 and Ksi= 0.5616).

myode[200, 0.201, 0.5616]

plt3 = 
      Plot[Evaluate[S[t] /. myode[200, 0.201, 0.5616]], {t, 0, 
          350}]

Show[plt3, plt2]