Phase-Plane Analysis of Biochemical Reactors
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Acknowledgement :
Numerical values were obtained from Process Dynamics, Modeling, Analysis and Simulation, B. Wayne Bequette, Prentice Hall, 1998.

Off[General::"spell1"]

Specific growth coefficient following form proposed by Monod in 1942.

(Subscript[μ, mod]) = (Subscript[K, mod]) x2/((Subscript[k, mod]) + x2)

(Subscript[K, mod]) = 0.53; (Subscript[k, mod]) = 0.12;

plt1 = 
      Plot[(Subscript[μ, mod]), {x2, 0, 4}, 
        PlotStyle -> RGBColor[1, 0, 0]]

Expression of the specific growth coefficient when there is substrate inhibition.

(Subscript[μ, SI]) = (Subscript[K, SI]) x2/((Subscript[k, SI]) + x2 + (Subscript[k, l]) x2^2)

(Subscript[K, SI]) = 0.53; (Subscript[k, SI]) = 0.12; (Subscript[k, l]) = 0.4545;

plt2 = 
      Plot[(Subscript[μ, SI]), {x2, 0, 4}, 
        PlotStyle -> RGBColor[0, 0, 1]]

Show[plt1, plt2]

Dilution rate is De=0.3; Yield is Y=0.4; x2f is the substrate concentration in the feed.

De = 0.3;

Y = 0.4;

x2f = 4.0;

We find two steady-state solutions for the Monod growth rate.

steady1 = 
      Solve[{((Subscript[μ, mod]) - De) x1 ==  0, 
          De (x2f - x2) - (Subscript[μ, mod]) x1/Y ==  0}, {x1, 
          x2}]

We find three steady-state solutions for the growth rate with substrate inhibition.

steady2 = 
      Solve[{((Subscript[μ, SI]) - De) x1 ==  0, 
          De (x2f - x2) - (Subscript[μ, SI]) x1/Y ==  0}, {x1, 
          x2}]

Phase-plane analysis when there is substrate inhibition in the growth rate. Red, Blue and Green dots in the figure are the three steady-states. The Red and Blue dots correspond to the stable steady-states. The Green dot is the saddle point (unstable steady-state). Only solution three (Blue dot) is non trivial. The Red dot corresponds to wash out (i.e. no biomass is left in the bioreactor).

(Subscript[μ, SI]) =. ;

(Subscript[μ, SI])[t_] := 
      (Subscript[K, SI]) x2[t]/((Subscript[k, SI]) + x2[t] + (Subscript[k, l]) x2[t]^2)

sol3 = 
      NDSolve[{D[
              x1[t], {t, 1}] ==  ((Subscript[μ, SI])[t] - De) x1[
                t], 
D[x2[t], {t, 1}] ==  
            De (x2f - x2[t]) - (Subscript[μ, SI])[t] x1[t]/
                  Y, 
x1[0] ==  1, 
          x2[0] ==  1}, {x1[t], x2[t]}, {t, 0, 30}]

plt5 = 
      Plot[{Evaluate[x1[t] /. sol3], Evaluate[x2[t] /. sol3]}, {t, 0, 
          30}, 
PlotStyle -> {RGBColor[1, 0, 0], 
            RGBColor[0, 0, 1]}]

sol4 = 
      NDSolve[{D[
              x1[t], {t, 1}] ==  ((Subscript[μ, SI])[t] - De) x1[
                t], 
D[x2[t], {t, 1}] ==  
            De (x2f - x2[t]) - (Subscript[μ, SI])[t] x1[t]/
                  Y, 
x1[0] ==  0.75, 
          x2[0] ==  2}, {x1[t], x2[t]}, {t, 0, 30}]

plt6 = 
      Plot[{Evaluate[x1[t] /. sol4], Evaluate[x2[t] /. sol4]}, {t, 0, 
          30}, 
PlotStyle -> {RGBColor[1, 0, 0], 
            RGBColor[0, 0, 1]}]

sys[x10_, 
        x20_] := {D[
            x1[t], {t, 1}] ==  ((Subscript[μ, SI])[t] - De) x1[
              t], 
D[x2[t], {t, 1}] ==  
          De (x2f - x2[t]) - (Subscript[μ, SI])[t] x1[t]/
                Y, 
x1[0] ==  x10, 
        x2[0] == x20}

sol[a_, b_] := 
      NDSolve[sys[x10, x20] /. x10 -> a /. x20 -> b, {x1[t], 
          x2[t]}, {t, 0, 30}]

For[i = 1, i < 4, i += 1, 
      plt[i] = ParametricPlot[
          Evaluate[{x1[t], x2[t]} /. sol[ 1.5, i]], {t, 0, 30}, 
          PlotRange -> All, 
          DisplayFunction -> Identity]]

For[i = 1, i < 4, i += 1, 
      PLT[i] = ParametricPlot[
          Evaluate[{x1[t], x2[t]} /. sol[0.1, i]], {t, 0, 30}, 
          PlotRange -> All, 
          DisplayFunction -> Identity]]

For[i = 1, i < 16, i += 1, 
      plot[i] = 
        ParametricPlot[
          Evaluate[{x1[t], x2[t]} /. sol[0.1 i, 0]], {t, 0, 30}, 
          PlotRange -> All, 
          DisplayFunction -> Identity]]

For[i = 1, i < 16, i += 1, 
      plotx[i] = 
        ParametricPlot[
          Evaluate[{x1[t], x2[t]} /. sol[0.1 i, 5]], {t, 0, 30}, 
          PlotRange -> All, 
          DisplayFunction -> Identity]]

tbl=Table[plt[i],{i,1,3}];

TBL=Table[PLT[i],{i,1,3}];

TbL=Table[plot[i],{i,1,15}];

TbLx = Table[plotx[i], {i, 1, 15}];

Show[Flatten[Join[TBL,tbl,TbL,TbLx]],Graphics[{{PointSize[0.03],RGBColor[1,0,0],Point[{0,4}]},{PointSize[0.03],RGBColor[0,1,0],Point[{0.9951,1.5123}]},{PointSize[0.03],RGBColor[0,0,1],Point[{1.5302,0.1745}]}}],PlotRange->{{0,2},{0,5}},DisplayFunction->$DisplayFunction]

Phase-plane analysis when the growth rate follows Monod equation. Red and Blue dots in the figure are the two steady-states. The Blue dot corresponds to the stable steady-state. The Red dot is the saddle point (unstable steady-state). Only solution two (Blue dot) is non trivial. The Red dot corresponds to wash out (i.e. no biomass is left in the bioreactor).

(Subscript[μ, mod]) =. ;

(Subscript[μ, mod])[t_] := (Subscript[K, mod]) x2[t]/((Subscript[k, mod]) + x2[t])

sol1 = 
      NDSolve[{D[
              x1[t], {t, 1}] ==  ((Subscript[μ, mod])[t] - De) x1[
                t], 
D[x2[t], {t, 1}] ==  
            De (x2f - x2[t]) - (Subscript[μ, mod])[t] x1[t]/
                  Y, 
x1[0] ==  1, 
          x2[0] ==  1}, {x1[t], x2[t]}, {t, 0, 30}]

plt3 = 
      Plot[{Evaluate[x1[t] /. sol1], Evaluate[x2[t] /. sol1]}, {t, 0, 
          30}, 
PlotStyle -> {RGBColor[1, 0, 0], 
            RGBColor[0, 0, 1]}]

sol2 = 
      NDSolve[{D[
              x1[t], {t, 1}] ==  ((Subscript[μ, mod])[t] - De) x1[
                t], 
D[x2[t], {t, 1}] ==  
            De (x2f - x2[t]) - (Subscript[μ, mod])[t] x1[t]/
                  Y, 
x1[0] ==  0.75, 
          x2[0] ==  2}, {x1[t], x2[t]}, {t, 0, 30}]

plt4 = 
      Plot[{Evaluate[x1[t] /. sol2], Evaluate[x2[t] /. sol2]}, {t, 0, 
          30}, 
PlotStyle -> {RGBColor[1, 0, 0], 
            RGBColor[0, 0, 1]}]

sys[x10_, 
        x20_] := {D[
            x1[t], {t, 1}] ==  ((Subscript[μ, mod])[t] - De) x1[
              t], 
D[x2[t], {t, 1}] ==  
          De (x2f - x2[t]) - (Subscript[μ, mod])[t] x1[t]/
                Y, 
x1[0] ==  x10, 
        x2[0] == x20}

sol[a_, b_] := 
      NDSolve[sys[x10, x20] /. x10 -> a /. x20 -> b, {x1[t], 
          x2[t]}, {t, 0, 30}]

For[i = 1, i < 16, i += 1, 
      plot[i] = 
        ParametricPlot[
          Evaluate[{x1[t], x2[t]} /. sol[0.1 i, 0]], {t, 0, 30}, 
          PlotRange -> All, 
          DisplayFunction -> Identity]]

For[i = 1, i < 16, i += 1, 
      plotx[i] = 
        ParametricPlot[
          Evaluate[{x1[t], x2[t]} /. sol[0.1 i, 5]], {t, 0, 30}, 
          PlotRange -> All, 
          DisplayFunction -> Identity]]

TbL=Table[plot[i],{i,1,10}];

TbLx = Table[plotx[i], {i, 1, 10}];

Show[Flatten[Join[TbL,TbLx]],Graphics[{{PointSize[0.03],RGBColor[1,0,0],Point[{0,4}]},{PointSize[0.03],RGBColor[0,0,1],Point[{1.53739,0.15652173}]}}],PlotRange->{{0,2},{0,5}},DisplayFunction->$DisplayFunction]

Stability analysis when the growth rate follows Monod equation. We look at the signs of the eigenvalues of the Jacobian Matrix.  Steady-state (0,4) is unstable because eigenvalues are of opposite sign indicating that it is a saddle point. Steady-state (1.53739,0.156522) is a stable steady-state because both eigenvalues are negative.

(Subscript[μ, mod]) = (Subscript[K, mod]) x2/((Subscript[k, mod]) + x2);

Eigen_steday1 = {{D[((Subscript[μ, mod]) - De) x1, {x1, 1}], 
              D[((Subscript[μ, mod]) - De) x1, {x2, 1}]}, {D[
                De (x2f - x2) - (Subscript[μ, mod]) x1/Y, {x1, 1}], 
              D[De (x2f - x2) - (Subscript[μ, mod]) x1/Y, {x2, 1}]}} /. 
          steady1[[1]] // Eigenvalues

Eigen_steday2 = {{D[((Subscript[μ, mod]) - De) x1, {x1, 1}], 
              D[((Subscript[μ, mod]) - De) x1, {x2, 1}]}, {D[
                De (x2f - x2) - (Subscript[μ, mod]) x1/Y, {x1, 1}], 
              D[De (x2f - x2) - (Subscript[μ, mod]) x1/Y, {x2, 1}]}} /. 
          steady1[[2]] // Eigenvalues

Stability analysis when the growth rate is subject to substrate inhibition. We look at the signs of the eigenvalues of the Jacobian Matrix.  Steady-states (0,4) and (1.53016,0.174593) are both stable because their eigenvalues are negative. Steady-state (0.995103,1.51224) is a saddle point because its eigenvalues have opposite signs.

(Subscript[μ, SI]) = (Subscript[K, SI]) x2/((Subscript[k, SI]) + x2 + (Subscript[k, l]) x2^2);

Eigen_steday1 = {{D[((Subscript[μ, SI]) - De) x1, {x1, 1}], 
              D[((Subscript[μ, SI]) - De) x1, {x2, 1}]}, {D[
                De (x2f - x2) - (Subscript[μ, SI]) x1/Y, {x1, 1}], 
              D[De (x2f - x2) - (Subscript[μ, SI]) x1/Y, {x2, 1}]}} /. 
          steady2[[3]] // Eigenvalues

Eigen_steday3 = {{D[((Subscript[μ, mod]) - De) x1, {x1, 1}], 
              D[((Subscript[μ, mod]) - De) x1, {x2, 1}]}, {D[
                De (x2f - x2) - (Subscript[μ, mod]) x1/Y, {x1, 1}], 
              D[De (x2f - x2) - (Subscript[μ, mod]) x1/Y, {x2, 1}]}} /. 
          steady2[[5]] // Eigenvalues

Eigen_steday2 = {{D[((Subscript[μ, mod]) - De) x1, {x1, 1}], 
              D[((Subscript[μ, mod]) - De) x1, {x2, 1}]}, {D[
                De (x2f - x2) - (Subscript[μ, mod]) x1/Y, {x1, 1}], 
              D[De (x2f - x2) - (Subscript[μ, mod]) x1/Y, {x2, 1}]}} /. 
          steady2[[4]] // Eigenvalues